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

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

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

Context Navigation

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

Last change on this file since 11250 was 11250, checked in by agn, 5 years ago
correctly specify dstokes for nn_osm_wave 0 & 1
Property svn:keywords set to `Id`
File size: 101.9 KB

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

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

Download in other formats: