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

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

source: NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90 @ 14316

Last change on this file since 14316 was 14316, checked in by smueller, 2 years ago

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