source: NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90 @ 14565

Last change on this file since 14565 was 14565, checked in by smueller, 7 months ago

Synchronisation of the OSMOSIS boundary layer scheme with the version developed in branch /NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0: transfer of changesets [14515,14519,14536] and completion of the transfer of changeset [14405] (ticket #2353)

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