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

Last change on this file since 14868 was 14868, checked in by smueller, 5 months ago

Removal of redundant halo exchanges in subroutine zdf_osm of module zdfosm (ticket #2353)

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