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

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

Bugfixes related to diagnostic output and various minor adjustments (ticket #2353)

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