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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfosm.F90 in NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90 @ 14404

Last change on this file since 14404 was 14404, checked in by agn, 3 years ago

Use alternate formulation for zsc_wth_pyc and zsc_ws_pyc

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