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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90 @ 9106

Last change on this file since 9106 was 9106, checked in by cetlod, 6 years ago
Bugfix on ZDF
File size: 88.5 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
	27	!! (14) Bouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out).
	28	!! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added.
	29	!! (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out)
	30	!! (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out)
	31	!!----------------------------------------------------------------------
[8946]	32
[8930]	33	!!----------------------------------------------------------------------
	34	!! 'key_zdfosm' OSMOSIS scheme
	35	!!----------------------------------------------------------------------
	36	!! zdf_osm : update momentum and tracer Kz from osm scheme
	37	!! zdf_osm_init : initialization, namelist read, and parameters control
	38	!! osm_rst : read (or initialize) and write osmosis restart fields
	39	!! tra_osm : compute and add to the T & S trend the non-local flux
	40	!! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD)
	41	!! dyn_osm : compute and add to u & v trensd the non-local flux
	42	!!----------------------------------------------------------------------
[8946]	43	USE oce ! ocean dynamics and active tracers
[8930]	44	! uses wn from previous time step (which is now wb) to calculate hbl
	45	USE dom_oce ! ocean space and time domain
	46	USE zdf_oce ! ocean vertical physics
	47	USE sbc_oce ! surface boundary condition: ocean
	48	USE sbcwave ! surface wave parameters
	49	USE phycst ! physical constants
	50	USE eosbn2 ! equation of state
	51	USE traqsr ! details of solar radiation absorption
	52	USE zdfddm ! double diffusion mixing (avs array)
	53	USE iom ! I/O library
	54	USE lib_mpp ! MPP library
	55	USE trd_oce ! ocean trends definition
	56	USE trdtra ! tracers trends
	57	!
	58	USE in_out_manager ! I/O manager
	59	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	60	USE prtctl ! Print control
	61	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
	62
	63	IMPLICIT NONE
	64	PRIVATE
	65
	66	PUBLIC zdf_osm ! routine called by step.F90
	67	PUBLIC zdf_osm_init ! routine called by nemogcm.F90
	68	PUBLIC osm_rst ! routine called by step.F90
	69	PUBLIC tra_osm ! routine called by step.F90
	70	PUBLIC trc_osm ! routine called by trcstp.F90
	71	PUBLIC dyn_osm ! routine called by 'step.F90'
	72
	73	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux
	74	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamv !: non-local v-momentum flux
	75	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamt !: non-local temperature flux (gamma/<ws>o)
	76	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghams !: non-local salinity flux (gamma/<ws>o)
	77	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt
	78	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth
[8946]	79	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbli !: intial boundary layer depth for stable blayer
	80	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift.
[8930]	81
	82	! !! Namelist namzdf_osm
	83	LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la
	84	REAL(wp) :: rn_osm_la ! Turbulent Langmuir number
	85	REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift
	86	REAL(wp) :: rn_osm_hbl0 = 10._wp ! Initial value of hbl for 1D runs
	87	INTEGER :: nn_ave ! = 0/1 flag for horizontal average on avt
	88	INTEGER :: nn_osm_wave = 0 ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
	89	LOGICAL :: ln_dia_osm ! Use namelist rn_osm_la
	90
	91
	92	LOGICAL :: ln_kpprimix = .true. ! Shear instability mixing
	93	REAL(wp) :: rn_riinfty = 0.7 ! local Richardson Number limit for shear instability
	94	REAL(wp) :: rn_difri = 0.005 ! maximum shear mixing at Rig = 0 (m2/s)
	95	LOGICAL :: ln_convmix = .true. ! Convective instability mixing
	96	REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s)
	97
	98	! !!! General constants
	99	REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number
	100	REAL(wp) :: pthird = 1._wp/3._wp ! 1/3
	101	REAL(wp) :: p2third = 2._wp/3._wp ! 2/3
	102
	103	INTEGER :: idebug = 236
	104	INTEGER :: jdebug = 228
	105	!!----------------------------------------------------------------------
	106	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
	107	!! $Id$
	108	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	109	!!----------------------------------------------------------------------
	110	CONTAINS
	111
	112	INTEGER FUNCTION zdf_osm_alloc()
	113	!!----------------------------------------------------------------------
	114	!! * FUNCTION zdf_osm_alloc *
	115	!!----------------------------------------------------------------------
	116	ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), &
	117	& hbl(jpi,jpj) , hbli(jpi,jpj) , dstokes(jpi, jpj) , &
	118	& etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
	119	IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
	120	IF( lk_mpp ) CALL mpp_sum ( zdf_osm_alloc )
	121	END FUNCTION zdf_osm_alloc
	122
[8946]	123
[8930]	124	SUBROUTINE zdf_osm( kt, p_avm, p_avt )
	125	!!----------------------------------------------------------------------
	126	!! * ROUTINE zdf_osm *
	127	!!
	128	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
	129	!! coefficients and non local mixing using the OSMOSIS scheme
	130	!!
	131	!! ** Method : The boundary layer depth hosm is diagnosed at tracer points
	132	!! from profiles of buoyancy, and shear, and the surface forcing.
	133	!! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
	134	!!
	135	!! Kx = hosm Wx(sigma) G(sigma)
	136	!!
	137	!! and the non local term ghamt = Cs / Ws(sigma) / hosm
	138	!! Below hosm the coefficients are the sum of mixing due to internal waves
	139	!! shear instability and double diffusion.
	140	!!
	141	!! -1- Compute the now interior vertical mixing coefficients at all depths.
	142	!! -2- Diagnose the boundary layer depth.
	143	!! -3- Compute the now boundary layer vertical mixing coefficients.
	144	!! -4- Compute the now vertical eddy vicosity and diffusivity.
	145	!! -5- Smoothing
	146	!!
	147	!! N.B. The computation is done from jk=2 to jpkm1
	148	!! Surface value of avt are set once a time to zero
	149	!! in routine zdf_osm_init.
	150	!!
	151	!! ** Action : update the non-local terms ghamts
	152	!! update avt (before vertical eddy coef.)
	153	!!
	154	!! References : Large W.G., Mc Williams J.C. and Doney S.C.
	155	!! Reviews of Geophysics, 32, 4, November 1994
	156	!! Comments in the code refer to this paper, particularly
	157	!! the equation number. (LMD94, here after)
	158	!!----------------------------------------------------------------------
[8946]	159	INTEGER , INTENT(in ) :: kt ! ocean time step
[8930]	160	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
	161	!!
	162	INTEGER :: ji, jj, jk ! dummy loop indices
	163	INTEGER :: ikbot, jkmax, jkm1, jkp2 !
	164
	165	REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube !
	166	REAL(wp) :: zalbet, zbeta, zthermal, zatt1 !
	167	REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
	168	REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm !
	169	REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density
	170	INTEGER :: jm ! dummy loop indices
	171	REAL(wp) :: zr1, zr2, zr3, zr4, zrhop ! Compression terms
	172	REAL(wp) :: zflag, zrn2, zdep21, zdep32, zdep43
	173	REAL(wp) :: zesh2, zri, zfri ! Interior richardson mixing
	174	REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
	175	REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages
	176	! Scales
	177	REAL(wp), DIMENSION(jpi,jpj) :: zrad0 ! Surface solar temperature flux (deg m/s)
	178	REAL(wp), DIMENSION(jpi,jpj) :: zradh ! Radiative flux at bl base (Buoyancy units)
	179	REAL(wp), DIMENSION(jpi,jpj) :: zradav ! Radiative flux, bl average (Buoyancy Units)
	180	REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity
	181	REAL(wp), DIMENSION(jpi,jpj) :: zwstrl ! Langmuir velocity scale
	182	REAL(wp), DIMENSION(jpi,jpj) :: zvstr ! Velocity scale that ends to zustar for large Langmuir number.
	183	REAL(wp), DIMENSION(jpi,jpj) :: zwstrc ! Convective velocity scale
	184	REAL(wp), DIMENSION(jpi,jpj) :: zuw0 ! Surface u-momentum flux
	185	REAL(wp), DIMENSION(jpi,jpj) :: zvw0 ! Surface v-momentum flux
	186	REAL(wp), DIMENSION(jpi,jpj) :: zwth0 ! Surface heat flux (Kinematic)
	187	REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux
	188	REAL(wp), DIMENSION(jpi,jpj) :: zwb0 ! Surface buoyancy flux
	189	REAL(wp), DIMENSION(jpi,jpj) :: zwthav ! Heat flux - bl average
	190	REAL(wp), DIMENSION(jpi,jpj) :: zwsav ! freshwater flux - bl average
	191	REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average
	192	REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux
	193	REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift
	194	REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number
	195	REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
	196	REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
	197	REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer
	198	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl
	199
	200	! mixed-layer variables
	201
	202	INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
	203	INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
	204
	205	REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
	206	REAL(wp) :: zugrad,zvgrad ! temporary variables for calculating pycnocline shear
	207
	208	REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid
	209	REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid
	210	REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid
	211	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
	212	REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer
	213	REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer
	214	REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer
	215	REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
	216	REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
	217	REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline
	218	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline
	219	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline
	220	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline
	221	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc ! u-shear across the pycnocline
	222	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc ! v-shear across the pycnocline
	223
	224	! Flux-gradient relationship variables
	225
	226	REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale.
	227
	228	REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity
	229	REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
	230	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.
	231	REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
	232
	233	! For calculating Ri#-dependent mixing
	234	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du ! u-shear^2
	235	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv ! v-shear^2
	236	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
	237
	238	! Temporary variables
	239	INTEGER :: inhml
	240	INTEGER :: i_lconv_alloc
	241	REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
	242	REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables
	243	REAL(wp) :: zthick, zz0, zz1 ! temporary variables
	244	REAL(wp) :: zvel_max, zhbl_s ! temporary variables
	245	REAL(wp) :: zfac ! temporary variable
	246	REAL(wp) :: zus_x, zus_y ! temporary Stokes drift
	247	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
	248	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
	249
	250	! For debugging
	251	INTEGER :: ikt
	252	!!--------------------------------------------------------------------
	253	!
	254	ALLOCATE( lconv(jpi,jpj), STAT= i_lconv_alloc )
	255	IF( i_lconv_alloc /= 0 ) CALL ctl_warn('zdf_osm: failed to allocate lconv')
	256
	257	ibld(:,:) = 0 ; imld(:,:) = 0
	258	zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp
	259	zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:) = 0._wp ; zuw0(:,:) = 0._wp
	260	zvw0(:,:) = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:) = 0._wp ; zwb0(:,:) = 0._wp
	261	zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:) = 0._wp ; zwb_ent(:,:) = 0._wp
	262	zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
	263	zhol(:,:) = 0._wp
	264	lconv(:,:) = .FALSE.
	265	! mixed layer
	266	! no initialization of zhbl or zhml (or zdh?)
	267	zhbl(:,:) = 1._wp ; zhml(:,:) = 1._wp ; zdh(:,:) = 1._wp ; zdhdt(:,:) = 0._wp
	268	zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp
	269	zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp
	270	zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp
	271	zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp
	272	zdt_ml(:,:) = 0._wp ; zds_ml(:,:) = 0._wp ; zdu_ml(:,:) = 0._wp ; zdv_ml(:,:) = 0._wp
	273	zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:) = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp
	274	zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp
	275	!
	276	zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
	277	zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
	278	!
	279	! Flux-Gradient arrays.
	280	zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp
	281	zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp
	282	zsc_wth_1(:,:) = 0._wp ; zsc_ws_1(:,:) = 0._wp ; zsc_uw_1(:,:) = 0._wp
	283	zsc_uw_2(:,:) = 0._wp ; zsc_vw_1(:,:) = 0._wp ; zsc_vw_2(:,:) = 0._wp
	284	zhbl_t(:,:) = 0._wp ; zdhdt(:,:) = 0._wp
	285
	286	zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
	287	ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp
	288
	289	! hbl = MAX(hbl,epsln)
	290	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	291	! Calculate boundary layer scales
	292	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	293
	294	! Assume two-band radiation model for depth of OSBL
	295	zz0 = rn_abs ! surface equi-partition in 2-bands
	296	zz1 = 1. - rn_abs
	297	DO jj = 2, jpjm1
	298	DO ji = 2, jpim1
	299	! Surface downward irradiance (so always +ve)
	300	zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp
	301	! Downwards irradiance at base of boundary layer
	302	zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
	303	! Downwards irradiance averaged over depth of the OSBL
	304	zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
	305	& + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
	306	END DO
	307	END DO
	308	! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
	309	DO jj = 2, jpjm1
	310	DO ji = 2, jpim1
	311	zthermal = rab_n(ji,jj,1,jp_tem)
	312	zbeta = rab_n(ji,jj,1,jp_sal)
	313	! Upwards surface Temperature flux for non-local term
	314	zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
	315	! Upwards surface salinity flux for non-local term
	316	zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
	317	! Non radiative upwards surface buoyancy flux
	318	zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj)
	319	! turbulent heat flux averaged over depth of OSBL
	320	zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
	321	! turbulent salinity flux averaged over depth of the OBSL
	322	zwsav(ji,jj) = 0.5 * zws0(ji,jj)
	323	! turbulent buoyancy flux averaged over the depth of the OBSBL
	324	zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj)
	325	! Surface upward velocity fluxes
	326	zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
	327	zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
	328	! Friction velocity (zustar), at T-point : LMD94 eq. 2
	329	zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
	330	zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
	331	zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
	332	END DO
	333	END DO
	334	! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
	335	SELECT CASE (nn_osm_wave)
	336	! Assume constant La#=0.3
	337	CASE(0)
	338	DO jj = 2, jpjm1
	339	DO ji = 2, jpim1
	340	zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
	341	zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
	342	zustke(ji,jj) = MAX ( SQRT( zus_xzus_x + zus_yzus_y), 1.0e-8 )
	343	! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
	344	END DO
	345	END DO
	346	! Assume Pierson-Moskovitz wind-wave spectrum
	347	CASE(1)
	348	DO jj = 2, jpjm1
	349	DO ji = 2, jpim1
	350	! Use wind speed wndm included in sbc_oce module
	351	zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
	352	dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav
	353	END DO
	354	END DO
	355	! Use ECMWF wave fields as output from SBCWAVE
	356	CASE(2)
	357	zfac = 2.0_wp * rpi / 16.0_wp
	358	DO jj = 2, jpjm1
	359	DO ji = 2, jpim1
	360	! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
	361	! The coefficient 0.8 gives La=0.3 in this situation.
	362	! It could represent the effects of the spread of wave directions
	363	! around the mean wind. The effect of this adjustment needs to be tested.
	364	zustke(ji,jj) = MAX (1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), zustar(ji,jj)/(0.45*0.45 ))
	365	dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)hsw(ji,jj) / ( MAX(zustke(ji,jj)wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
	366	END DO
	367	END DO
	368	END SELECT
	369
	370	! Langmuir velocity scale (zwstrl), La # (zla)
	371	! mixed scale (zvstr), convective velocity scale (zwstrc)
	372	DO jj = 2, jpjm1
	373	DO ji = 2, jpim1
	374	! Langmuir velocity scale (zwstrl), at T-point
	375	zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
	376	! Modify zwstrl to allow for small and large values of dstokes/hbl.
	377	! Intended as a possible test. Doesn't affect LES results for entrainment,
	378	! but hasn't been shown to be correct as dstokes/h becomes large or small.
	379	zwstrl(ji,jj) = zwstrl(ji,jj) * &
	380	& (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))*pthird &
	381	& ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) ))
	382	! define La this way so effects of Stokes penetration depth on velocity scale are included
	383	zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3
	384	! Velocity scale that tends to zustar for large Langmuir numbers
	385	zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + &
	386	& ( 1.0 - EXP( -0.5 * zla(ji,jj)*2 ) ) zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
	387
	388	! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
	389	! Note zustke and zwstrl are not amended.
	390	IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45
	391	!
	392	! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
	393	IF ( zwbav(ji,jj) > 0.0) THEN
	394	zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
	395	zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
	396	lconv(ji,jj) = .TRUE.
	397	ELSE
	398	zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln )
	399	lconv(ji,jj) = .FALSE.
	400	ENDIF
	401	END DO
	402	END DO
	403
	404	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	405	! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
	406	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	407	! BL must be always 2 levels deep.
	408	hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) )
	409	ibld(:,:) = 3
	410	DO jk = 4, jpkm1
	411	DO jj = 2, jpjm1
	412	DO ji = 2, jpim1
	413	IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
	414	ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
	415	ENDIF
	416	END DO
	417	END DO
	418	END DO
	419
	420	DO jj = 2, jpjm1 ! Vertical slab
	421	DO ji = 2, jpim1
	422	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
	423	zbeta = rab_n(ji,jj,1,jp_sal)
	424	zt = 0._wp
	425	zs = 0._wp
	426	zu = 0._wp
	427	zv = 0._wp
	428	! average over depth of boundary layer
	429	zthick=0._wp
	430	DO jm = 2, ibld(ji,jj)
	431	zthick=zthick+e3t_n(ji,jj,jm)
	432	zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
	433	zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
	434	zu = zu + e3t_n(ji,jj,jm) &
	435	& * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
	436	& / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
	437	zv = zv + e3t_n(ji,jj,jm) &
	438	& * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
	439	& / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
	440	END DO
	441	zt_bl(ji,jj) = zt / zthick
	442	zs_bl(ji,jj) = zs / zthick
	443	zu_bl(ji,jj) = zu / zthick
	444	zv_bl(ji,jj) = zv / zthick
	445	zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
	446	zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
	447	zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
	448	& / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
	449	zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
	450	& / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
	451	zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
	452	IF ( lconv(ji,jj) ) THEN ! Convective
	453	zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) &
	454	& + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) )
	455
	456	zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
	457	& ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
	458	! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment.
	459	! zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) &
	460	! & + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)*2 ) zustar(ji,jj)**3/hbl(ji,jj) )
	461
	462	! zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / &
	463	! & ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
	464	zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) )
	465	ELSE ! Stable
	466	zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) &
	467	& + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) &
	468	& - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) &
	469	& * zwstrl(ji,jj)**3 / hbli(ji,jj)
	470	zzdhdt = zzdhdt + zwbav(ji,jj)
	471	IF ( zzdhdt < 0._wp ) THEN
	472	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
	473	zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj)
	474	ELSE
	475	zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) + MAX( zdb_bl(ji,jj),0.0 )
	476	ENDIF
	477	zzdhdt = 2.0 * zzdhdt / zpert
	478	ENDIF
	479	zdhdt(ji,jj) = zzdhdt
	480	END DO
	481	END DO
	482
	483	! Calculate averages over depth of boundary layer
	484	imld = ibld ! use imld to hold previous blayer index
	485	ibld(:,:) = 3
	486
	487	zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it
	488	zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:))
	489	zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom
	490
	491	DO jk = 4, jpkm1
	492	DO jj = 2, jpjm1
	493	DO ji = 2, jpim1
	494	IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
	495	ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
	496	ENDIF
	497	END DO
	498	END DO
	499	END DO
	500
	501	!
	502	! Step through model levels taking account of buoyancy change to determine the effect on dhdt
	503	!
	504	DO jj = 2, jpjm1
	505	DO ji = 2, jpim1
	506	IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
	507	!
	508	! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
	509	!
	510	zhbl_s = hbl(ji,jj)
	511	jm = imld(ji,jj)
	512	zthermal = rab_n(ji,jj,1,jp_tem)
	513	zbeta = rab_n(ji,jj,1,jp_sal)
	514	IF ( lconv(ji,jj) ) THEN
	515	!unstable
	516	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) / &
	517	& ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
	518
	519	DO jk = imld(ji,jj), ibld(ji,jj)
	520	zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
	521	& - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max
	522
	523	zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) )
	524	zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
	525
	526	IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
	527	END DO
	528	hbl(ji,jj) = zhbl_s
	529	ibld(ji,jj) = jm
	530	hbli(ji,jj) = hbl(ji,jj)
	531	ELSE
	532	! stable
	533	DO jk = imld(ji,jj), ibld(ji,jj)
	534	zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + &
	535	& 2.0 * zwstrl(ji,jj)**2 / zhbl_s
	536
	537	zhbl_s = zhbl_s + (&
	538	& 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) + &
	539	& ( ( 0.32 / 3.0 ) * EXP( -2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) - &
	540	& ( 0.32 / 3.0 - 0.0485 ) * EXP( -12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) * zwstrl(ji,jj)*3 / hbli(ji,jj) ) / zdb &
	541	& e3w_n(ji,jj,jk) / zdhdt(ji,jj) ! ALMG to investigate whether need to include wn here
	542
	543	zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
	544	IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
	545	END DO
	546	hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) )
	547	ibld(ji,jj) = MAX(jm, 3 )
	548	IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
	549	ENDIF ! IF ( lconv )
	550	ELSE
	551	! change zero or one model level.
	552	hbl(ji,jj) = zhbl_t(ji,jj)
	553	IF ( lconv(ji,jj) ) THEN
	554	hbli(ji,jj) = hbl(ji,jj)
	555	ELSE
	556	hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) )
	557	IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
	558	ENDIF
	559	ENDIF
	560	zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
	561	END DO
	562	END DO
	563	dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms.
	564
	565	! Recalculate averages over boundary layer after depth updated
	566	! Consider later combining this into the loop above and looking for columns
	567	! where the index for base of the boundary layer have changed
	568	DO jj = 2, jpjm1 ! Vertical slab
	569	DO ji = 2, jpim1
	570	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
	571	zbeta = rab_n(ji,jj,1,jp_sal)
	572	zt = 0._wp
	573	zs = 0._wp
	574	zu = 0._wp
	575	zv = 0._wp
	576	! average over depth of boundary layer
	577	zthick=0._wp
	578	DO jm = 2, ibld(ji,jj)
	579	zthick=zthick+e3t_n(ji,jj,jm)
	580	zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
	581	zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
	582	zu = zu + e3t_n(ji,jj,jm) &
	583	& * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
	584	& / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
	585	zv = zv + e3t_n(ji,jj,jm) &
	586	& * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
	587	& / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
	588	END DO
	589	zt_bl(ji,jj) = zt / zthick
	590	zs_bl(ji,jj) = zs / zthick
	591	zu_bl(ji,jj) = zu / zthick
	592	zv_bl(ji,jj) = zv / zthick
	593	zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
	594	zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
	595	zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
	596	& / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
	597	zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
	598	& / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
	599	zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
	600	zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
	601	IF ( lconv(ji,jj) ) THEN
	602	IF ( zdb_bl(ji,jj) > 0._wp )THEN
	603	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability
	604	zari = 4.5 * ( zvstr(ji,jj)**2 ) &
	605	& / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01
	606	ELSE ! unstable
	607	zari = 4.5 * ( zwstrc(ji,jj)**2 ) &
	608	& / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01
	609	ENDIF
	610	IF ( zari > 0.2 ) THEN ! This test checks for weakly stratified pycnocline
	611	zari = 0.2
	612	zwb_ent(ji,jj) = 0._wp
	613	ENDIF
	614	inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
	615	imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
	616	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
	617	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
	618	ELSE ! IF (zdb_bl)
	619	imld(ji,jj) = ibld(ji,jj) - 1
	620	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
	621	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
	622	ENDIF
	623	ELSE ! IF (lconv)
	624	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
	625	! boundary layer deepening
	626	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
	627	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
	628	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
	629	& / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 )
	630	inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
	631	imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
	632	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
	633	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
	634	ELSE
	635	imld(ji,jj) = ibld(ji,jj) - 1
	636	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
	637	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
	638	ENDIF ! IF (zdb_bl > 0.0)
	639	ELSE ! IF(dhdt >= 0)
	640	! boundary layer collapsing.
	641	imld(ji,jj) = ibld(ji,jj)
	642	zhml(ji,jj) = zhbl(ji,jj)
	643	zdh(ji,jj) = 0._wp
	644	ENDIF ! IF (dhdt >= 0)
	645	ENDIF ! IF (lconv)
	646	END DO
	647	END DO
	648
	649	! Average over the depth of the mixed layer in the convective boundary layer
	650	! Also calculate entrainment fluxes for temperature and salinity
	651	DO jj = 2, jpjm1 ! Vertical slab
	652	DO ji = 2, jpim1
	653	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
	654	zbeta = rab_n(ji,jj,1,jp_sal)
	655	IF ( lconv(ji,jj) ) THEN
	656	zt = 0._wp
	657	zs = 0._wp
	658	zu = 0._wp
	659	zv = 0._wp
	660	! average over depth of boundary layer
	661	zthick=0._wp
	662	DO jm = 2, imld(ji,jj)
	663	zthick=zthick+e3t_n(ji,jj,jm)
	664	zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
	665	zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
	666	zu = zu + e3t_n(ji,jj,jm) &
	667	& * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
	668	& / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
	669	zv = zv + e3t_n(ji,jj,jm) &
	670	& * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
	671	& / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
	672	END DO
	673	zt_ml(ji,jj) = zt / zthick
	674	zs_ml(ji,jj) = zs / zthick
	675	zu_ml(ji,jj) = zu / zthick
	676	zv_ml(ji,jj) = zv / zthick
	677	zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
	678	zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
	679	zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
	680	& / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
	681	zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
	682	& / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
	683	zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
	684	ELSE
	685	! stable, if entraining calulate average below interface layer.
	686	IF ( zdhdt(ji,jj) >= 0._wp ) THEN
	687	zt = 0._wp
	688	zs = 0._wp
	689	zu = 0._wp
	690	zv = 0._wp
	691	! average over depth of boundary layer
	692	zthick=0._wp
	693	DO jm = 2, imld(ji,jj)
	694	zthick=zthick+e3t_n(ji,jj,jm)
	695	zt = zt + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
	696	zs = zs + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
	697	zu = zu + e3t_n(ji,jj,jm) &
	698	& * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
	699	& / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
	700	zv = zv + e3t_n(ji,jj,jm) &
	701	& * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
	702	& / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
	703	END DO
	704	zt_ml(ji,jj) = zt / zthick
	705	zs_ml(ji,jj) = zs / zthick
	706	zu_ml(ji,jj) = zu / zthick
	707	zv_ml(ji,jj) = zv / zthick
	708	zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem)
	709	zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
	710	zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
	711	& / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
	712	zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
	713	& / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
	714	zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
	715	ENDIF
	716	ENDIF
	717	END DO
	718	END DO
	719	!
	720	! rotate mean currents and changes onto wind align co-ordinates
	721	!
	722
	723	DO jj = 2, jpjm1
	724	DO ji = 2, jpim1
	725	ztemp = zu_ml(ji,jj)
	726	zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj)
	727	zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
	728	ztemp = zdu_ml(ji,jj)
	729	zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj)
	730	zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
	731	!
	732	ztemp = zu_bl(ji,jj)
	733	zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj)
	734	zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
	735	ztemp = zdu_bl(ji,jj)
	736	zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj)
	737	zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
	738	END DO
	739	END DO
	740
	741	zuw_bse = 0._wp
	742	zvw_bse = 0._wp
	743	DO jj = 2, jpjm1
	744	DO ji = 2, jpim1
	745
	746	IF ( lconv(ji,jj) ) THEN
	747	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
	748	zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
	749	zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
	750	ENDIF
	751	ELSE
	752	zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )
	753	zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )
	754	ENDIF
	755	END DO
	756	END DO
	757
	758	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	759	! Pycnocline gradients for scalars and velocity
	760	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	761
	762	DO jj = 2, jpjm1
	763	DO ji = 2, jpim1
	764	!
	765	IF ( lconv (ji,jj) ) THEN
	766	! Unstable conditions
	767	IF( zdb_bl(ji,jj) > 0._wp ) THEN
	768	! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero
	769	ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) )
	770	zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) )
	771	zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) )
	772	DO jk = 2 , ibld(ji,jj)
	773	znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
	774	zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	775	zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	776	zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	777	END DO
	778	ENDIF
	779	ELSE
	780	! stable conditions
	781	! if pycnocline profile only defined when depth steady of increasing.
	782	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady.
	783	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
	784	IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline
	785	ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
	786	zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
	787	zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
	788	DO jk = 2, ibld(ji,jj)
	789	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	790	zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
	791	zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
	792	zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
	793	END DO
	794	ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
	795	ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
	796	zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
	797	zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
	798	DO jk = 2, ibld(ji,jj)
	799	znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
	800	zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	801	zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	802	zdsdz_pyc(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	803	END DO
	804	ENDIF ! IF (zhol >=0.5)
	805	ENDIF ! IF (zdb_bl> 0.)
	806	ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero
	807	ENDIF ! IF (lconv)
	808	!
	809	END DO
	810	END DO
	811	!
	812	DO jj = 2, jpjm1
	813	DO ji = 2, jpim1
	814	!
	815	IF ( lconv (ji,jj) ) THEN
	816	! Unstable conditions
	817	zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / &
	818	& (( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0)
	819	zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / &
	820	& ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
	821	DO jk = 2 , ibld(ji,jj)-1
	822	znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
	823	zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	824	zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	825	END DO
	826	ELSE
	827	! stable conditions
	828	zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
	829	zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
	830	DO jk = 2, ibld(ji,jj)
	831	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	832	IF ( znd < 1.0 ) THEN
	833	zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
	834	ELSE
	835	zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
	836	ENDIF
	837	zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
	838	END DO
	839	ENDIF
	840	!
	841	END DO
	842	END DO
	843	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	844	! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
	845	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	846
	847	! WHERE ( lconv )
	848	! zdifml_sc = zhml * ( zwstrl*3 + 0.5 zwstrc3 )pthird
	849	! zvisml_sc = zdifml_sc
	850	! zdifpyc_sc = 0.165 * ( zwstrl3 + zwstrc3 )*pthird ( zhbl - zhml )
	851	! zvispyc_sc = 0.142 * ( zwstrl*3 + 0.5 zwstrc3 )pthird * ( zhbl - zhml )
	852	! zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third
	853	! zbeta_v_sc = 1.0 - 2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml
	854	! ELSEWHERE
	855	! zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )
	856	! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )
	857	! ENDWHERE
	858	DO jj = 2, jpjm1
	859	DO ji = 2, jpim1
	860	IF ( lconv(ji,jj) ) THEN
	861	zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
	862	zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
	863	zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zdh(ji,jj)
	864	zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zdh(ji,jj)
	865	zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third
	866	zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj)
	867	ELSE
	868	zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
	869	zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
	870	END IF
	871	END DO
	872	END DO
	873	!
	874	DO jj = 2, jpjm1
	875	DO ji = 2, jpim1
	876	IF ( lconv(ji,jj) ) THEN
	877	DO jk = 2, imld(ji,jj) ! mixed layer diffusivity
	878	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
	879	!
	880	zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
	881	!
	882	zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) * ( 1.0 - 0.5 * zznd_ml**2 )
	883	END DO
	884	! pycnocline - if present linear profile
	885	IF ( zdh(ji,jj) > 0._wp ) THEN
	886	DO jk = imld(ji,jj)+1 , ibld(ji,jj)
	887	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
	888	!
	889	zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj)*( 1.0 + zznd_pyc )
	890	!
	891	zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
	892	END DO
	893	ENDIF
	894	! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out
	895	zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj))
	896	! could be taken out, take account of entrainment represents as a diffusivity
	897	! should remove w from here, represents entrainment
	898	ELSE
	899	! stable conditions
	900	DO jk = 2, ibld(ji,jj)
	901	zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	902	zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
	903	zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
	904	END DO
	905	ENDIF ! end if ( lconv )
	906	!
	907	END DO ! end of ji loop
	908	END DO ! end of jj loop
	909
	910	!
	911	! calculate non-gradient components of the flux-gradient relationships
	912	!
	913	! Stokes term in scalar flux, flux-gradient relationship
	914	WHERE ( lconv )
	915	zsc_wth_1 = zwstrl*3 zwth0 / ( zvstr*3 + 0.5 zwstrc**3 + epsln)
	916	!
	917	zsc_ws_1 = zwstrl*3 zws0 / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
	918	ELSEWHERE
	919	zsc_wth_1 = 2.0 * zwthav
	920	!
	921	zsc_ws_1 = 2.0 * zwsav
	922	ENDWHERE
	923
	924
	925	DO jj = 2, jpjm1
	926	DO ji = 2, jpim1
	927	IF ( lconv(ji,jj) ) THEN
	928	DO jk = 2, imld(ji,jj)
	929	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	930	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)
	931	!
	932	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)
	933	END DO ! end jk loop
	934	ELSE ! else for if (lconv)
	935	! Stable conditions
	936	DO jk = 2, ibld(ji,jj)
	937	zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	938	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
	939	!
	940	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
	941	END DO
	942	ENDIF ! endif for check on lconv
	943
	944	END DO ! end of ji loop
	945	END DO ! end of jj loop
	946
	947
	948	! 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)
	949	WHERE ( lconv )
	950	zsc_uw_1 = ( zwstrl*3 + 0.5 zwstrc3 )pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) )
	951	zsc_uw_2 = ( zwstrl*3 + 0.5 zwstrc3 )pthird * zustke / ( zla**(8.0/3.0) + epsln )
	952	zsc_vw_1 = ff_t * zhml * zustke*3 zla(8.0/3.0) / ( ( zvstr3 + 0.5 * zwstrc3 )(2.0/3.0) + epsln )
	953	ELSEWHERE
	954	zsc_uw_1 = zustar**2
	955	zsc_vw_1 = ff_t * zhbl * zustke*3 zla(8.0/3.0) / (zvstr2 + epsln)
	956	ENDWHERE
	957
	958	DO jj = 2, jpjm1
	959	DO ji = 2, jpim1
	960	IF ( lconv(ji,jj) ) THEN
	961	DO jk = 2, imld(ji,jj)
	962	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	963	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) + 0.00125 * EXP ( -zznd_d ) * zsc_uw_2(ji,jj) ) * ( 1.0 - EXP ( -2.0 * zznd_d ) )
	964	!
	965	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
	966	END DO ! end jk loop
	967	ELSE
	968	! Stable conditions
	969	DO jk = 2, ibld(ji,jj) ! corrected to ibld
	970	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	971	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
	972	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
	973	END DO ! end jk loop
	974	ENDIF
	975	END DO ! ji loop
	976	END DO ! jj loo
	977
	978	! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
	979
	980	WHERE ( lconv )
	981	zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
	982	zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
	983	ELSEWHERE
	984	zsc_wth_1 = 0._wp
	985	zsc_ws_1 = 0._wp
	986	ENDWHERE
	987
	988	DO jj = 2, jpjm1
	989	DO ji = 2, jpim1
	990	IF (lconv(ji,jj) ) THEN
	991	DO jk = 2, imld(ji,jj)
	992	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
	993	! calculate turbulent length scale
	994	zl_c = 0.9 * ( 1.0 - EXP ( -7.0 * ( zznd_ml - zznd_ml*3 / 3.0 ) ) ) ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) )
	995	zl_l = 2.0 * ( 1.0 - EXP ( -2.0 * ( zznd_ml - zznd_ml*3 / 3.0 ) ) ) (1.0 - EXP ( -5.0 * ( 1.0 -zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
	996	zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
	997	! non-gradient buoyancy terms
	998	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 )
	999	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 )
	1000	END DO
	1001	ELSE
	1002	DO jk = 2, ibld(ji,jj)
	1003	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
	1004	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj)
	1005	END DO
	1006	ENDIF
	1007	END DO ! ji loop
	1008	END DO ! jj loop
	1009
	1010
	1011	WHERE ( lconv )
	1012	zsc_uw_1 = -zwb0 * zustar*2 zhml / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
	1013	zsc_uw_2 = zwb0 * zustke * zhml / ( zvstr*3 + 0.5 zwstrc3 + epsln )(2.0/3.0)
	1014	zsc_vw_1 = 0._wp
	1015	ELSEWHERE
	1016	zsc_uw_1 = 0._wp
	1017	zsc_vw_1 = 0._wp
	1018	ENDWHERE
	1019
	1020	DO jj = 2, jpjm1
	1021	DO ji = 2, jpim1
	1022	IF ( lconv(ji,jj) ) THEN
	1023	DO jk = 2 , imld(ji,jj)
	1024	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	1025	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) * ( 1.0 - EXP ( -0.5 * zznd_d) ) * zsc_uw_2(ji,jj) )
	1026	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
	1027	END DO ! jk loop
	1028	ELSE
	1029	! stable conditions
	1030	DO jk = 2, ibld(ji,jj)
	1031	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
	1032	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
	1033	END DO
	1034	ENDIF
	1035	END DO ! ji loop
	1036	END DO ! jj loop
	1037
	1038	! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
	1039
	1040	WHERE ( lconv )
	1041	zsc_wth_1 = zwth0
	1042	zsc_ws_1 = zws0
	1043	ELSEWHERE
	1044	zsc_wth_1 = 2.0 * zwthav
	1045	zsc_ws_1 = zws0
	1046	ENDWHERE
	1047
	1048	DO jj = 2, jpjm1
	1049	DO ji = 2, jpim1
	1050	IF ( lconv(ji,jj) ) THEN
	1051	DO jk = 2, imld(ji,jj)
	1052	zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
	1053	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml*4 ) - EXP ( -6.0 zznd_ml ) ) ) &
	1054	& * ( 1.0 -EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
	1055	!
	1056	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml*4 ) - EXP ( -6.0 zznd_ml ) ) ) &
	1057	& * ( 1.0 -EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
	1058	END DO
	1059	ELSE
	1060	DO jk = 2, ibld(ji,jj)
	1061	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	1062	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	1063	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
	1064	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
	1065	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
	1066	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
	1067	END DO
	1068	ENDIF
	1069	ENDDO ! ji loop
	1070	END DO ! jj loop
	1071
	1072
	1073	WHERE ( lconv )
	1074	zsc_uw_1 = zustar**2
	1075	zsc_vw_1 = ff_t * zustke * zhml
	1076	ELSEWHERE
	1077	zsc_uw_1 = zustar**2
	1078	zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
	1079	zsc_vw_1 = ff_t * zustke * zhbl
	1080	zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )*2 ) zsc_vw_1
	1081	ENDWHERE
	1082
	1083	DO jj = 2, jpjm1
	1084	DO ji = 2, jpim1
	1085	IF ( lconv(ji,jj) ) THEN
	1086	DO jk = 2, imld(ji,jj)
	1087	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
	1088	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	1089	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
	1090	& + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml*4 ) - EXP ( -8.0 zznd_ml ) ) * zsc_uw_1(ji,jj)
	1091	!
	1092	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
	1093	& + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
	1094	END DO
	1095	ELSE
	1096	DO jk = 2, ibld(ji,jj)
	1097	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	1098	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
	1099	IF ( zznd_d <= 2.0 ) THEN
	1100	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
	1101	&* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
	1102	!
	1103	ELSE
	1104	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
	1105	& + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
	1106	!
	1107	ENDIF
	1108
	1109	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
	1110	& + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d*2 ) zsc_vw_1(ji,jj)
	1111	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
	1112	& + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
	1113	END DO
	1114	ENDIF
	1115	END DO
	1116	END DO
	1117	!
	1118	! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
	1119
	1120	DO jj = 2, jpjm1
	1121	DO ji = 2, jpim1
	1122	IF ( lconv(ji,jj) ) THEN
	1123	DO jk = 2, ibld(ji,jj)
	1124	znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
	1125	IF ( znd >= 0.0 ) THEN
	1126	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
	1127	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
	1128	ELSE
	1129	ghamu(ji,jj,jk) = 0._wp
	1130	ghamv(ji,jj,jk) = 0._wp
	1131	ENDIF
	1132	END DO
	1133	ELSE
	1134	DO jk = 2, ibld(ji,jj)
	1135	znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
	1136	IF ( znd >= 0.0 ) THEN
	1137	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
	1138	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
	1139	ELSE
	1140	ghamu(ji,jj,jk) = 0._wp
	1141	ghamv(ji,jj,jk) = 0._wp
	1142	ENDIF
	1143	END DO
	1144	ENDIF
	1145	END DO
	1146	END DO
	1147
	1148	! pynocline contributions
	1149	! Temporary fix to avoid instabilities when zdb_bl becomes very very small
	1150	zsc_uw_1 = 0._wp ! 50.0 * zla*(8.0/3.0) zustar*2 zhbl / ( zdb_bl + epsln )
	1151	DO jj = 2, jpjm1
	1152	DO ji = 2, jpim1
	1153	DO jk= 2, ibld(ji,jj)
	1154	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	1155	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
	1156	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
	1157	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
	1158	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )*(7.0/4.0) zdbdz_pyc(ji,jj,jk)
	1159	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
	1160	END DO
	1161	END DO
	1162	END DO
	1163
	1164	! Entrainment contribution.
	1165
	1166	DO jj=2, jpjm1
	1167	DO ji = 2, jpim1
	1168	IF ( lconv(ji,jj) ) THEN
	1169	DO jk = 1, imld(ji,jj) - 1
	1170	znd=gdepw_n(ji,jj,jk) / zhml(ji,jj)
	1171	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
	1172	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
	1173	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd
	1174	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd
	1175	END DO
	1176	DO jk = imld(ji,jj), ibld(ji,jj)
	1177	znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
	1178	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
	1179	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
	1180	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd )
	1181	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd )
	1182	END DO
	1183	ENDIF
	1184	ghamt(ji,jj,ibld(ji,jj)) = 0._wp
	1185	ghams(ji,jj,ibld(ji,jj)) = 0._wp
	1186	ghamu(ji,jj,ibld(ji,jj)) = 0._wp
	1187	ghamv(ji,jj,ibld(ji,jj)) = 0._wp
	1188	END DO ! ji loop
	1189	END DO ! jj loop
	1190
	1191
	1192	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	1193	! Need to put in code for contributions that are applied explicitly to
	1194	! the prognostic variables
	1195	! 1. Entrainment flux
	1196	!
	1197	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	1198
	1199
	1200
	1201	! rotate non-gradient velocity terms back to model reference frame
	1202
	1203	DO jj = 2, jpjm1
	1204	DO ji = 2, jpim1
	1205	DO jk = 2, ibld(ji,jj)
	1206	ztemp = ghamu(ji,jj,jk)
	1207	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
	1208	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
	1209	END DO
	1210	END DO
	1211	END DO
	1212
	1213	IF(ln_dia_osm) THEN
	1214	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
	1215	END IF
	1216
	1217	! KPP-style Ri# mixing
	1218	IF( ln_kpprimix) THEN
	1219	DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form)
	1220	DO jj = 1, jpjm1
	1221	DO ji = 1, jpim1 ! vector opt.
	1222	z3du(ji,jj,jk) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) &
	1223	& * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) &
	1224	& / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
	1225	z3dv(ji,jj,jk) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) &
	1226	& * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) &
	1227	& / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
	1228	END DO
	1229	END DO
	1230	END DO
	1231	!
	1232	DO jk = 2, jpkm1
	1233	DO jj = 2, jpjm1
	1234	DO ji = 2, jpim1 ! vector opt.
	1235	! ! shear prod. at w-point weightened by mask
	1236	zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
	1237	& + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
	1238	! ! local Richardson number
	1239	zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
	1240	zfri = MIN( zri / rn_riinfty , 1.0_wp )
	1241	zfri = ( 1.0_wp - zfri * zfri )
	1242	zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk)
	1243	END DO
	1244	END DO
	1245	END DO
	1246
	1247	DO jj = 2, jpjm1
	1248	DO ji = 2, jpim1
	1249	DO jk = ibld(ji,jj) + 1, jpkm1
	1250	zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
	1251	zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
	1252	END DO
	1253	END DO
	1254	END DO
	1255
	1256	END IF ! ln_kpprimix = .true.
	1257
	1258	! KPP-style set diffusivity large if unstable below BL
	1259	IF( ln_convmix) THEN
	1260	DO jj = 2, jpjm1
	1261	DO ji = 2, jpim1
	1262	DO jk = ibld(ji,jj) + 1, jpkm1
	1263	IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
	1264	END DO
	1265	END DO
	1266	END DO
	1267	END IF ! ln_convmix = .true.
	1268
	1269	! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
	1270	CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
	1271
	1272	! GN 25/8: need to change tmask --> wmask
	1273
	1274	DO jk = 2, jpkm1
	1275	DO jj = 2, jpjm1
	1276	DO ji = 2, jpim1
	1277	p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
	1278	p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
	1279	END DO
	1280	END DO
	1281	END DO
	1282	! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids
[9106]	1283	CALL lbc_lnk_multi( p_avt, 'W', 1. , p_avm, 'W', 1., &
[9104]	1284	& ghamu, 'W', 1. , ghamv, 'W', 1. )
[8930]	1285	DO jk = 2, jpkm1
	1286	DO jj = 2, jpjm1
	1287	DO ji = 2, jpim1
	1288	ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
	1289	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
	1290
	1291	ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
	1292	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
	1293
	1294	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk)
	1295	ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk)
	1296	END DO
	1297	END DO
	1298	END DO
	1299	! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged)
[9104]	1300	! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged)
[9106]	1301	CALL lbc_lnk_multi( ghamt, 'W', 1. , ghams, 'W', 1., &
[9104]	1302	& ghamu, 'U', 1. , ghamv, 'V', 1. )
[8930]	1303
	1304	IF(ln_dia_osm) THEN
	1305	SELECT CASE (nn_osm_wave)
	1306	! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1).
	1307	CASE(0:1)
	1308	IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)zustkezcos_wind ) ! x surface Stokes drift
	1309	IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)zustkezsin_wind ) ! y surface Stokes drift
	1310	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2zustke )
	1311	! Stokes drift read in from sbcwave (=2).
	1312	CASE(2)
	1313	IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd ) ! x surface Stokes drift
	1314	IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd ) ! y surface Stokes drift
	1315	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2 &
	1316	& SQRT(ut0sd2 + vt0sd2 ) )
	1317	END SELECT
	1318	IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL>
	1319	IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL>
	1320	IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL>
	1321	IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL>
	1322	IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! <Tw_0>
	1323	IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0>
	1324	IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth
	1325	IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli ) ! Initial boundary-layer depth
	1326	IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth
	1327	IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points
	1328	IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc ) ! convective velocity scale
	1329	IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale
	1330	IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale
	1331	IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.rau0tmask(:,:,1)zustar*3 ) ! BL depth internal to zdf_osm routine
	1332	IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.rau0tmask(:,:,1)zustar2zustke )
	1333	IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine
	1334	IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine
	1335	IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! ML depth internal to zdf_osm routine
	1336	IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine
	1337	IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! ML depth internal to zdf_osm routine
	1338	IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! ML depth internal to zdf_osm routine
	1339	IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML
[8946]	1340	END IF
	1341	! Lateral boundary conditions on p_avt (sign unchanged)
	1342	CALL lbc_lnk( p_avt(:,:,:), 'W', 1. )
	1343	!
[8930]	1344	END SUBROUTINE zdf_osm
	1345
[8946]	1346
[8930]	1347	SUBROUTINE zdf_osm_init
	1348	!!----------------------------------------------------------------------
	1349	!! * ROUTINE zdf_osm_init *
	1350	!!
	1351	!! ** Purpose : Initialization of the vertical eddy diffivity and
	1352	!! viscosity when using a osm turbulent closure scheme
	1353	!!
	1354	!! ** Method : Read the namosm namelist and check the parameters
	1355	!! called at the first timestep (nit000)
	1356	!!
	1357	!! ** input : Namlist namosm
	1358	!!----------------------------------------------------------------------
	1359	INTEGER :: ios ! local integer
	1360	INTEGER :: ji, jj, jk ! dummy loop indices
	1361	!!
	1362	NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
	1363	& ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
	1364	& ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv
	1365	!!----------------------------------------------------------------------
	1366	!
	1367	REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model
	1368	READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
[9104]	1369	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist', lwp )
[8930]	1370
	1371	REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
	1372	READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
[9104]	1373	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist', lwp )
[8930]	1374	IF(lwm) WRITE ( numond, namzdf_osm )
	1375
	1376	IF(lwp) THEN ! Control print
	1377	WRITE(numout,*)
	1378	WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
	1379	WRITE(numout,*) '~~~~~~~~~~~~'
	1380	WRITE(numout,*) ' Namelist namzdf_osm : set tke mixing parameters'
	1381	WRITE(numout,*) ' Use namelist rn_osm_la ln_use_osm_la = ', ln_use_osm_la
	1382	WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la
	1383	WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0
	1384	WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes
	1385	WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave
	1386	WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave
	1387	SELECT CASE (nn_osm_wave)
	1388	CASE(0)
	1389	WRITE(numout,*) ' calculated assuming constant La#=0.3'
	1390	CASE(1)
	1391	WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves'
	1392	CASE(2)
	1393	WRITE(numout,*) ' calculated from ECMWF wave fields'
	1394	END SELECT
	1395	WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm
	1396	WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix
	1397	WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
	1398	WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri
	1399	WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix
	1400	WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv
	1401	ENDIF
	1402
	1403	! ! allocate zdfosm arrays
	1404	IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
	1405
	1406	call osm_rst( nit000, 'READ' ) !* read or initialize hbl
	1407
	1408	IF( ln_zdfddm) THEN
	1409	IF(lwp) THEN
	1410	WRITE(numout,*)
	1411	WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
	1412	WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm '
	1413	ENDIF
	1414	ENDIF
	1415
	1416
	1417	!set constants not in namelist
	1418	!-----------------------------
	1419
	1420	IF(lwp) THEN
	1421	WRITE(numout,*)
	1422	ENDIF
	1423
	1424	IF (nn_osm_wave == 0) THEN
	1425	dstokes(:,:) = rn_osm_dstokes
	1426	END IF
	1427
	1428	! Horizontal average : initialization of weighting arrays
	1429	! -------------------
	1430
	1431	SELECT CASE ( nn_ave )
	1432
	1433	CASE ( 0 ) ! no horizontal average
	1434	IF(lwp) WRITE(numout,*) ' no horizontal average on avt'
	1435	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
	1436	! weighting mean arrays etmean
	1437	! ( 1 1 )
	1438	! avt = 1/4 ( 1 1 )
	1439	!
	1440	etmean(:,:,:) = 0.e0
	1441
	1442	DO jk = 1, jpkm1
	1443	DO jj = 2, jpjm1
	1444	DO ji = 2, jpim1 ! vector opt.
	1445	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
	1446	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
	1447	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
	1448	END DO
	1449	END DO
	1450	END DO
	1451
	1452	CASE ( 1 ) ! horizontal average
	1453	IF(lwp) WRITE(numout,*) ' horizontal average on avt'
	1454	! weighting mean arrays etmean
	1455	! ( 1/2 1 1/2 )
	1456	! avt = 1/8 ( 1 2 1 )
	1457	! ( 1/2 1 1/2 )
	1458	etmean(:,:,:) = 0.e0
	1459
	1460	DO jk = 1, jpkm1
	1461	DO jj = 2, jpjm1
	1462	DO ji = 2, jpim1 ! vector opt.
	1463	etmean(ji,jj,jk) = tmask(ji, jj,jk) &
	1464	& / MAX( 1., 2.* tmask(ji,jj,jk) &
	1465	& +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) &
	1466	& +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
	1467	& +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) &
	1468	& +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) )
	1469	END DO
	1470	END DO
	1471	END DO
	1472
	1473	CASE DEFAULT
	1474	WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave
	1475	CALL ctl_stop( ctmp1 )
	1476
	1477	END SELECT
	1478
	1479	! Initialization of vertical eddy coef. to the background value
	1480	! -------------------------------------------------------------
	1481	DO jk = 1, jpk
	1482	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
	1483	END DO
	1484
	1485	! zero the surface flux for non local term and osm mixed layer depth
	1486	! ------------------------------------------------------------------
	1487	ghamt(:,:,:) = 0.
	1488	ghams(:,:,:) = 0.
	1489	ghamu(:,:,:) = 0.
	1490	ghamv(:,:,:) = 0.
	1491	!
	1492	END SUBROUTINE zdf_osm_init
	1493
[8946]	1494
[8930]	1495	SUBROUTINE osm_rst( kt, cdrw )
	1496	!!---------------------------------------------------------------------
	1497	!! * ROUTINE osm_rst *
	1498	!!
	1499	!! ** Purpose : Read or write BL fields in restart file
	1500	!!
	1501	!! ** Method : use of IOM library. If the restart does not contain
	1502	!! required fields, they are recomputed from stratification
	1503	!!----------------------------------------------------------------------
	1504
	1505	INTEGER, INTENT(in) :: kt
	1506	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	1507
	1508	INTEGER :: id1, id2 ! iom enquiry index
	1509	INTEGER :: ji, jj, jk ! dummy loop indices
	1510	INTEGER :: iiki, ikt ! local integer
	1511	REAL(wp) :: zhbf ! tempory scalars
	1512	REAL(wp) :: zN2_c ! local scalar
	1513	REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth
	1514	INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top)
	1515	!!----------------------------------------------------------------------
	1516	!
	1517	!!-----------------------------------------------------------------------------
	1518	! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
	1519	!!-----------------------------------------------------------------------------
	1520	IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
	1521	id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. )
	1522	IF( id1 > 0 ) THEN ! 'wn' exists; read
	1523	CALL iom_get( numror, jpdom_autoglo, 'wn', wn )
	1524	WRITE(numout,*) ' ===>>>> : wn read from restart file'
	1525	ELSE
	1526	wn(:,:,:) = 0._wp
	1527	WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
	1528	END IF
	1529	id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. )
	1530	id2 = iom_varid( numror, 'hbli' , ldstop = .FALSE. )
	1531	IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return
	1532	CALL iom_get( numror, jpdom_autoglo, 'hbl', hbl )
	1533	CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli )
	1534	WRITE(numout,*) ' ===>>>> : hbl & hbli read from restart file'
	1535	RETURN
	1536	ELSE ! 'hbl' & 'hbli' not in restart file, recalculate
	1537	WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
	1538	END IF
	1539	END IF
	1540
	1541	!!-----------------------------------------------------------------------------
	1542	! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
	1543	!!-----------------------------------------------------------------------------
	1544	IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return
	1545	IF(lwp) WRITE(numout,*) '---- osm-rst ----'
	1546	CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn )
	1547	CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl )
	1548	CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli )
	1549	RETURN
	1550	END IF
	1551
	1552	!!-----------------------------------------------------------------------------
	1553	! Getting hbl, no restart file with hbl, so calculate from surface stratification
	1554	!!-----------------------------------------------------------------------------
	1555	IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
	1556	ALLOCATE( imld_rst(jpi,jpj) )
	1557	! w-level of the mixing and mixed layers
	1558	CALL eos_rab( tsn, rab_n )
	1559	CALL bn2(tsn, rab_n, rn2)
	1560	imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point
	1561	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
	1562	zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
	1563	!
	1564	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
	1565	DO jk = 1, jpkm1
	1566	DO jj = 1, jpj ! Mixed layer level: w-level
	1567	DO ji = 1, jpi
	1568	ikt = mbkt(ji,jj)
	1569	hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
	1570	IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
	1571	END DO
	1572	END DO
	1573	END DO
	1574	!
	1575	DO jj = 1, jpj
	1576	DO ji = 1, jpi
	1577	iiki = imld_rst(ji,jj)
	1578	hbl (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth
	1579	END DO
	1580	END DO
	1581	hbl = MAX(hbl,epsln)
	1582	hbli(:,:) = hbl(:,:)
	1583	DEALLOCATE( imld_rst )
	1584	WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
	1585	END SUBROUTINE osm_rst
	1586
[8946]	1587
[8930]	1588	SUBROUTINE tra_osm( kt )
	1589	!!----------------------------------------------------------------------
	1590	!! * ROUTINE tra_osm *
	1591	!!
	1592	!! ** Purpose : compute and add to the tracer trend the non-local tracer flux
	1593	!!
	1594	!! ** Method : ???
	1595	!!----------------------------------------------------------------------
	1596	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace
	1597	!!----------------------------------------------------------------------
	1598	INTEGER, INTENT(in) :: kt
	1599	INTEGER :: ji, jj, jk
	1600	!
	1601	IF( kt == nit000 ) THEN
	1602	IF(lwp) WRITE(numout,*)
	1603	IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
	1604	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	1605	ENDIF
	1606
	1607	IF( l_trdtra ) THEN !* Save ta and sa trends
	1608	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
	1609	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)
	1610	ENDIF
	1611
	1612	! add non-local temperature and salinity flux
	1613	DO jk = 1, jpkm1
	1614	DO jj = 2, jpjm1
	1615	DO ji = 2, jpim1
	1616	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &
	1617	& - ( ghamt(ji,jj,jk ) &
	1618	& - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
	1619	tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) &
	1620	& - ( ghams(ji,jj,jk ) &
	1621	& - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
	1622	END DO
	1623	END DO
	1624	END DO
	1625
	1626
	1627	! save the non-local tracer flux trends for diagnostic
	1628	IF( l_trdtra ) THEN
	1629	ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
	1630	ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
	1631	!!bug gm jpttdzdf ==> jpttosm
	1632	CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
	1633	CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
	1634	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
	1635	ENDIF
	1636
	1637	IF(ln_ctl) THEN
	1638	CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm - Ta: ', mask1=tmask, &
	1639	& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
	1640	ENDIF
	1641	!
	1642	END SUBROUTINE tra_osm
	1643
[8946]	1644
[8930]	1645	SUBROUTINE trc_osm( kt ) ! Dummy routine
	1646	!!----------------------------------------------------------------------
	1647	!! * ROUTINE trc_osm *
	1648	!!
	1649	!! ** Purpose : compute and add to the passive tracer trend the non-local
	1650	!! passive tracer flux
	1651	!!
	1652	!!
	1653	!! ** Method : ???
	1654	!!----------------------------------------------------------------------
[8946]	1655	!
[8930]	1656	!!----------------------------------------------------------------------
	1657	INTEGER, INTENT(in) :: kt
	1658	WRITE(,) 'trc_osm: Not written yet', kt
	1659	END SUBROUTINE trc_osm
	1660
[8946]	1661
[8930]	1662	SUBROUTINE dyn_osm( kt )
	1663	!!----------------------------------------------------------------------
	1664	!! * ROUTINE dyn_osm *
	1665	!!
	1666	!! ** Purpose : compute and add to the velocity trend the non-local flux
	1667	!! copied/modified from tra_osm
	1668	!!
	1669	!! ** Method : ???
	1670	!!----------------------------------------------------------------------
[8946]	1671	INTEGER, INTENT(in) :: kt !
	1672	!
	1673	INTEGER :: ji, jj, jk ! dummy loop indices
[8930]	1674	!!----------------------------------------------------------------------
	1675	!
	1676	IF( kt == nit000 ) THEN
	1677	IF(lwp) WRITE(numout,*)
	1678	IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
	1679	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	1680	ENDIF
	1681	!code saving tracer trends removed, replace with trdmxl_oce
	1682
[8946]	1683	DO jk = 1, jpkm1 ! add non-local u and v fluxes
[8930]	1684	DO jj = 2, jpjm1
	1685	DO ji = 2, jpim1
	1686	ua(ji,jj,jk) = ua(ji,jj,jk) &
	1687	& - ( ghamu(ji,jj,jk ) &
	1688	& - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
	1689	va(ji,jj,jk) = va(ji,jj,jk) &
	1690	& - ( ghamv(ji,jj,jk ) &
	1691	& - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
	1692	END DO
	1693	END DO
	1694	END DO
[9089]	1695	!
[8930]	1696	! code for saving tracer trends removed
	1697	!
	1698	END SUBROUTINE dyn_osm
	1699
[8946]	1700	!!======================================================================
[8930]	1701	END MODULE zdfosm

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

Download in other formats: