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/UKMO/NEMO_4.0_OSMOSIS/src/OCE/ZDF – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/src/OCE/ZDF/zdfosm.F90 @ 11145

Last change on this file since 11145 was 11145, checked in by cguiavarch, 5 years ago
Merge changes from /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser:r11139-11143
File size: 101.8 KB

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

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

Download in other formats: