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/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfosm.F90 @ 14051

Last change on this file since 14051 was 14051, checked in by techene, 4 years ago

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