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

Last change on this file since 14889 was 14889, checked in by smueller, 8 months ago

Enabling of the extended halo option (nn_hls=2) in the implementation of the OSMOSIS boundary layer scheme (ticket #2353)

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