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

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

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

Last change on this file since 14824 was 14824, checked in by smueller, 3 years ago

Removal of optional debug output (ticket #2353)
  • Property svn:keywords set to Id
File size: 206.1 KB
Line 
1MODULE zdfosm
2   !!======================================================================
3   !!                       ***  MODULE  zdfosm  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the OSMOSIS
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !!  History : NEMO 4.0  ! A. Grant, G. Nurser
8   !! 15/03/2017  Changed calculation of pycnocline thickness in unstable conditions and stable conditions AG
9   !! 15/03/2017  Calculation of pycnocline gradients for stable conditions changed. Pycnocline gradients now depend on stability of the OSBL. A.G
10   !! 06/06/2017  (1) Checks on sign of buoyancy jump in calculation of  OSBL depth.  A.G.
11   !!             (2) Removed variable zbrad0, zbradh and zbradav since they are not used.
12   !!             (3) Approximate treatment for shear turbulence.
13   !!                        Minimum values for zustar and zustke.
14   !!                        Add velocity scale, zvstr, that tends to zustar for large Langmuir numbers.
15   !!                        Limit maximum value for Langmuir number.
16   !!                        Use zvstr in definition of stability parameter zhol.
17   !!             (4) Modified parametrization of entrainment flux, changing original coefficient 0.0485 for Langmuir contribution to 0.135 * zla
18   !!             (5) For stable boundary layer add factor that depends on length of timestep to 'slow' collapse and growth. Make sure buoyancy jump not negative.
19   !!             (6) For unstable conditions when growth is over multiple levels, limit change to maximum of one level per cycle through loop.
20   !!             (7) Change lower limits for loops that calculate OSBL averages from 1 to 2. Large gradients between levels 1 and 2 can cause problems.
21   !!             (8) Change upper limits from ibld-1 to ibld.
22   !!             (9) Calculation of pycnocline thickness in unstable conditions. Check added to ensure that buoyancy jump is positive before calculating Ri.
23   !!            (10) Thickness of interface layer at base of the stable OSBL set by Richardson number. Gives continuity in transition from unstable OSBL.
24   !!            (11) Checks that buoyancy jump is poitive when calculating pycnocline profiles.
25   !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity.
26   !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information
27   !!            (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence.
28   !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added.
29   !!            (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out)
30   !!            (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out)
31   !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made,
32   !!                  (a) Pycnocline temperature and salinity profies changed for unstable layers
33   !!                  (b) The stable OSBL depth parametrization changed.
34   !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code.
35   !! 23/05/19   (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1
36   !!----------------------------------------------------------------------
37
38   !!----------------------------------------------------------------------
39   !!   'ln_zdfosm'                                             OSMOSIS scheme
40   !!----------------------------------------------------------------------
41   !!   zdf_osm        : update momentum and tracer Kz from osm scheme
42   !!      zdf_osm_vertical_average             : compute vertical averages over boundary layers
43   !!      zdf_osm_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   USE timing, ONLY : timing_start, timing_stop   ! Timer
84   !
85   IMPLICIT NONE
86   PRIVATE
87   !
88   ! Public subroutines
89   PUBLIC zdf_osm        ! Routine called by step.F90
90   PUBLIC zdf_osm_init   ! Routine called by nemogcm.F90
91   PUBLIC osm_rst        ! Routine called by step.F90
92   PUBLIC tra_osm        ! Routine called by step.F90
93   PUBLIC trc_osm        ! Routine called by trcstp.F90
94   PUBLIC dyn_osm        ! Routine called by step.F90
95   !
96   ! Public variables
97   LOGICAL,  PUBLIC                                      ::   ln_osm_mle   !: Flag to activate the Mixed Layer Eddy (MLE)
98   !                                                                       !     parameterisation, needed by tra_mle_init in
99   !                                                                       !     tramle.F90
100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu        !: Non-local u-momentum flux
101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv        !: Non-local v-momentum flux
102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamt        !: Non-local temperature flux (gamma/<ws>o)
103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghams        !: Non-local salinity flux (gamma/<ws>o)
104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl          !: Boundary layer depth
105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml          !: ML depth
106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle         !: Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization
107   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle     !: Zonal buoyancy gradient in ML
108   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle     !: Meridional buoyancy gradient in ML
109   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof     !: Level of base of MLE layer
110   !
111   INTERFACE zdf_osm_velocity_rotation
112      !!---------------------------------------------------------------------
113      !!              ***  INTERFACE zdf_velocity_rotation  ***
114      !!---------------------------------------------------------------------
115      MODULE PROCEDURE zdf_osm_velocity_rotation_2d
116      MODULE PROCEDURE zdf_osm_velocity_rotation_3d
117   END INTERFACE
118   !
119   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   ! Averaging operator for avt
120   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! Depth of pycnocline
121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_ft    ! Inverse of the modified Coriolis parameter at t-pts
122   !
123   ! Layer indices
124   INTEGER,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nbld   ! Level of boundary layer base
125   INTEGER,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmld   ! Level of mixed-layer depth (pycnocline top)
126   !
127   ! Layer type
128   INTEGER,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   n_ddh   ! Type of shear layer
129   !                                                        !    n_ddh=0: active shear layer
130   !                                                        !    n_ddh=1: shear layer not active
131   !                                                        !    n_ddh=2: shear production low
132   !
133   ! Layer flags
134   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_conv    ! Unstable/stable bl
135   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_shear   ! Shear layers
136   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_coup    ! Coupling to bottom
137   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_pyc     ! OSBL pycnocline present
138   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_flux    ! Surface flux extends below OSBL into MLE layer
139   LOGICAL,  ALLOCATABLE, SAVE, DIMENSION(:,:) ::   l_mle     ! MLE layer increases in hickness.
140   !
141   ! Scales
142   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swth0       ! Surface heat flux (Kinematic)
143   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sws0        ! Surface freshwater flux
144   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swb0        ! Surface buoyancy flux
145   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   suw0        ! Surface u-momentum flux
146   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sustar      ! Friction velocity
147   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scos_wind   ! Cos angle of surface stress
148   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssin_wind   ! Sin angle of surface stress
149   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swthav      ! Heat flux - bl average
150   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swsav       ! Freshwater flux - bl average
151   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swbav       ! Buoyancy flux - bl average
152   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sustke      ! Surface Stokes drift
153   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dstokes     ! Penetration depth of the Stokes drift
154   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swstrl      ! Langmuir velocity scale
155   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   swstrc      ! Convective velocity scale
156   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sla         ! Trubulent Langmuir number
157   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   svstr       ! Velocity scale that tends to sustar for large Langmuir number
158   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shol        ! Stability parameter for boundary layer
159   !
160   ! Layer averages: BL
161   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_t_bl   ! Temperature average
162   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_s_bl   ! Salinity average
163   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_u_bl   ! Velocity average (u)
164   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_v_bl   ! Velocity average (v)
165   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_bl   ! Buoyancy average
166   !
167   ! Difference between layer average and parameter at the base of the layer: BL
168   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dt_bl   ! Temperature difference
169   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_ds_bl   ! Salinity difference
170   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_du_bl   ! Velocity difference (u)
171   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dv_bl   ! Velocity difference (v)
172   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_db_bl   ! Buoyancy difference
173   !
174   ! Layer averages: ML
175   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_t_ml   ! Temperature average
176   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_s_ml   ! Salinity average
177   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_u_ml   ! Velocity average (u)
178   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_v_ml   ! Velocity average (v)
179   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_ml   ! Buoyancy average
180   !
181   ! Difference between layer average and parameter at the base of the layer: ML
182   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dt_ml   ! Temperature difference
183   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_ds_ml   ! Salinity difference
184   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_du_ml   ! Velocity difference (u)
185   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_dv_ml   ! Velocity difference (v)
186   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_db_ml   ! Buoyancy difference
187   !
188   ! Layer averages: MLE
189   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_t_mle  ! Temperature average
190   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_s_mle  ! Salinity average
191   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_u_mle  ! Velocity average (u)
192   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_v_mle  ! Velocity average (v)
193   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   av_b_mle  ! Buoyancy average
194   !
195   ! Diagnostic output
196   REAL(WP), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   osmdia2d   ! Auxiliary array for diagnostic output
197   REAL(WP), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   osmdia3d   ! Auxiliary array for diagnostic output
198   !
199   !            ** Namelist  namzdf_osm  **
200   LOGICAL  ::   ln_use_osm_la                      ! Use namelist rn_osm_la
201   REAL(wp) ::   rn_osm_la                          ! Turbulent Langmuir number
202   REAL(wp) ::   rn_osm_dstokes                     ! Depth scale of Stokes drift
203   REAL(wp) ::   rn_zdfosm_adjust_sd   = 1.0_wp     ! Factor to reduce Stokes drift by
204   REAL(wp) ::   rn_osm_hblfrac        = 0.1_wp     ! For nn_osm_wave = 3/4 specify fraction in top of hbl
205   LOGICAL  ::   ln_zdfosm_ice_shelter              ! Flag to activate ice sheltering
206   REAL(wp) ::   rn_osm_hbl0           = 10.0_wp    ! Initial value of hbl for 1D runs
207   INTEGER  ::   nn_ave                             ! = 0/1 flag for horizontal average on avt
208   INTEGER  ::   nn_osm_wave = 0                    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into
209   !                                                !    sbcwave
210   INTEGER  ::   nn_osm_SD_reduce                   ! = 0/1/2 flag for getting effective stokes drift from surface value
211   LOGICAL  ::   ln_dia_osm                         ! Use namelist  rn_osm_la
212   LOGICAL  ::   ln_dia_pyc_scl        = .FALSE.    ! Output of pycnocline scalar-gradient profiles
213   LOGICAL  ::   ln_dia_pyc_shr        = .FALSE.    ! Output of pycnocline velocity-shear  profiles
214   LOGICAL  ::   ln_kpprimix           = .TRUE.     ! Shear instability mixing
215   REAL(wp) ::   rn_riinfty            = 0.7_wp     ! Local Richardson Number limit for shear instability
216   REAL(wp) ::   rn_difri              = 0.005_wp   ! Maximum shear mixing at Rig = 0    (m2/s)
217   LOGICAL  ::   ln_convmix            = .TRUE.     ! Convective instability mixing
218   REAL(wp) ::   rn_difconv            = 1.0_wp     ! Diffusivity when unstable below BL  (m2/s)
219   !
220   ! OSMOSIS mixed layer eddy parametrization constants
221   INTEGER  ::   nn_osm_mle          ! = 0/1 flag for horizontal average on avt
222   REAL(wp) ::   rn_osm_mle_ce       ! MLE coefficient
223   !   ! Parameters used in nn_osm_mle = 0 case
224   REAL(wp) ::   rn_osm_mle_lf       ! Typical scale of mixed layer front
225   REAL(wp) ::   rn_osm_mle_time     ! Time scale for mixing momentum across the mixed layer
226   !   ! Parameters used in nn_osm_mle = 1 case
227   REAL(wp) ::   rn_osm_mle_lat      ! Reference latitude for a 5 km scale of ML front
228   LOGICAL  ::   ln_osm_hmle_limit   ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
229   REAL(wp) ::   rn_osm_hmle_limit   ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
230   REAL(wp) ::   rn_osm_mle_rho_c    ! Density criterion for definition of MLD used by FK
231   REAL(wp) ::   rb_c                ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld
232   REAL(wp) ::   rc_f                ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case
233   REAL(wp) ::   rn_osm_mle_thresh   ! Threshold buoyancy for deepening of MLE layer below OSBL base
234   REAL(wp) ::   rn_osm_bl_thresh    ! Threshold buoyancy for deepening of OSBL base
235   REAL(wp) ::   rn_osm_mle_tau      ! Adjustment timescale for MLE
236   !
237   !             ** General constants  **
238   REAL(wp) ::   epsln     = 1.0e-20_wp      ! A small positive number to ensure no div by zero
239   REAL(wp) ::   depth_tol = 1.0e-6_wp       ! A small-ish positive number to give a hbl slightly shallower than gdepw
240   REAL(wp) ::   pthird    = 1.0_wp/3.0_wp   ! 1/3
241   REAL(wp) ::   p2third   = 2.0_wp/3.0_wp   ! 2/3
242   !
243   !! * Substitutions
244#  include "do_loop_substitute.h90"
245#  include "domzgr_substitute.h90"
246   !!----------------------------------------------------------------------
247   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
248   !! $Id$
249   !! Software governed by the CeCILL license (see ./LICENSE)
250   !!----------------------------------------------------------------------
251CONTAINS
252
253   INTEGER FUNCTION zdf_osm_alloc()
254      !!----------------------------------------------------------------------
255      !!                 ***  FUNCTION zdf_osm_alloc  ***
256      !!----------------------------------------------------------------------
257      INTEGER ::   ierr
258      !
259      zdf_osm_alloc = 0
260      !
261      ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), hbl(jpi,jpj), hml(jpi,jpj),   &
262         &      hmle(jpi,jpj),      dbdx_mle(jpi,jpj),  dbdy_mle(jpi,jpj),  mld_prof(jpi,jpj),  STAT=ierr )
263      zdf_osm_alloc = zdf_osm_alloc + ierr
264      !
265      ALLOCATE( etmean(A2D(0),jpk), dh(jpi,jpj), r1_ft(A2D(0)), STAT=ierr )
266      zdf_osm_alloc = zdf_osm_alloc + ierr
267      !
268      ALLOCATE( nbld(jpi,jpj), nmld(A2D(0)), STAT=ierr )
269      zdf_osm_alloc = zdf_osm_alloc + ierr
270      !
271      ALLOCATE( n_ddh(A2D(0)), STAT=ierr )
272      zdf_osm_alloc = zdf_osm_alloc + ierr
273      !
274      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 )
275      zdf_osm_alloc = zdf_osm_alloc + ierr
276      !
277      ALLOCATE( swth0(A2D(0)),     sws0(A2D(0)),   swb0(A2D(0)),  suw0(A2D(0)),  sustar(A2D(0)), scos_wind(A2D(0)),   &
278         &      ssin_wind(A2D(0)), swthav(A2D(0)), swsav(A2D(0)), swbav(A2D(0)), sustke(A2D(0)), dstokes(A2D(0)),     &
279         &      swstrl(A2D(0)),    swstrc(A2D(0)), sla(A2D(0)),   svstr(A2D(0)), shol(A2D(0)),   STAT=ierr )
280      zdf_osm_alloc = zdf_osm_alloc + ierr
281      !
282      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)
283      zdf_osm_alloc = zdf_osm_alloc + ierr
284      !
285      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)
286      zdf_osm_alloc = zdf_osm_alloc + ierr
287      !
288      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)
289      zdf_osm_alloc = zdf_osm_alloc + ierr
290      !
291      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)
292      zdf_osm_alloc = zdf_osm_alloc + ierr
293      !
294      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)
295      zdf_osm_alloc = zdf_osm_alloc + ierr
296      !
297      IF ( ln_dia_osm ) THEN
298         ALLOCATE( osmdia2d(jpi,jpj), osmdia3d(jpi,jpj,jpk), STAT=ierr )
299         zdf_osm_alloc = zdf_osm_alloc + ierr
300      END IF
301      !
302      CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
303      IF( zdf_osm_alloc /= 0 ) CALL ctl_warn( 'zdf_osm_alloc: failed to allocate zdf_osm arrays' )
304      !
305   END FUNCTION zdf_osm_alloc
306
307   SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm,   &
308      &                p_avt )
309      !!----------------------------------------------------------------------
310      !!                   ***  ROUTINE zdf_osm  ***
311      !!
312      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
313      !!      coefficients and non local mixing using the OSMOSIS scheme
314      !!
315      !! ** Method :   The boundary layer depth hosm is diagnosed at tracer points
316      !!      from profiles of buoyancy, and shear, and the surface forcing.
317      !!      Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
318      !!
319      !!                      Kx =  hosm  Wx(sigma) G(sigma)
320      !!
321      !!             and the non local term ghamt = Cs / Ws(sigma) / hosm
322      !!      Below hosm  the coefficients are the sum of mixing due to internal waves
323      !!      shear instability and double diffusion.
324      !!
325      !!      -1- Compute the now interior vertical mixing coefficients at all depths.
326      !!      -2- Diagnose the boundary layer depth.
327      !!      -3- Compute the now boundary layer vertical mixing coefficients.
328      !!      -4- Compute the now vertical eddy vicosity and diffusivity.
329      !!      -5- Smoothing
330      !!
331      !!        N.B. The computation is done from jk=2 to jpkm1
332      !!             Surface value of avt are set once a time to zero
333      !!             in routine zdf_osm_init.
334      !!
335      !! ** Action  :   update the non-local terms ghamts
336      !!                update avt (before vertical eddy coef.)
337      !!
338      !! References : Large W.G., Mc Williams J.C. and Doney S.C.
339      !!         Reviews of Geophysics, 32, 4, November 1994
340      !!         Comments in the code refer to this paper, particularly
341      !!         the equation number. (LMD94, here after)
342      !!----------------------------------------------------------------------
343      INTEGER                   , INTENT(in   ) ::  kt               ! Ocean time step
344      INTEGER                   , INTENT(in   ) ::  Kbb, Kmm, Krhs   ! Ocean time level indices
345      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt     ! Momentum and tracer Kz (w-points)
346      !
347      ! Local variables
348      INTEGER ::   ji, jj, jk, jl, jm, jkflt   ! Dummy loop indices
349      !
350      REAL(wp) ::   zthermal, zbeta
351      REAL(wp) ::   zesh2, zri, zfri   ! Interior Richardson mixing
352      !
353      ! Scales
354      REAL(wp), DIMENSION(A2D(0))  ::   zrad0       ! Surface solar temperature flux (deg m/s)
355      REAL(wp), DIMENSION(A2D(0))  ::   zradh       ! Radiative flux at bl base (Buoyancy units)
356      REAL(wp)                     ::   zradav      ! Radiative flux, bl average (Buoyancy Units)
357      REAL(wp)                     ::   zvw0        ! Surface v-momentum flux
358      REAL(wp), DIMENSION(A2D(0))  ::   zwb0tot     ! Total surface buoyancy flux including insolation
359      REAL(wp), DIMENSION(A2D(0))  ::   zwb_ent     ! Buoyancy entrainment flux
360      REAL(wp), DIMENSION(A2D(0))  ::   zwb_min
361      REAL(wp), DIMENSION(A2D(0))  ::   zwb_fk_b    ! MLE buoyancy flux averaged over OSBL
362      REAL(wp), DIMENSION(A2D(0))  ::   zwb_fk      ! Max MLE buoyancy flux
363      REAL(wp), DIMENSION(A2D(0))  ::   zdiff_mle   ! Extra MLE vertical diff
364      REAL(wp), DIMENSION(A2D(0))  ::   zvel_mle    ! Velocity scale for dhdt with stable ML and FK
365      !
366      ! mixed-layer variables
367      INTEGER,  DIMENSION(A2D(0)) ::   jp_ext   ! Offset for external level
368      !
369      REAL(wp), DIMENSION(A2D(0)) ::   zhbl   ! BL depth - grid
370      REAL(wp), DIMENSION(A2D(0)) ::   zhml   ! ML depth - grid
371      !
372      REAL(wp), DIMENSION(A2D(0))  ::   zhmle   ! MLE depth - grid
373      REAL(wp), DIMENSION(jpi,jpj) ::   zmld   ! ML depth on grid
374      !
375      REAL(wp), DIMENSION(A2D(0))  ::   zdh   ! Pycnocline depth - grid
376      REAL(wp), DIMENSION(A2D(0))  ::   zdhdt   ! BL depth tendency
377      REAL(wp), DIMENSION(A2D(0))  ::   zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext   ! External temperature/salinity and buoyancy gradients
378      REAL(wp), DIMENSION(jpi,jpj) ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization
379      !
380      REAL(wp), DIMENSION(A2D(0)) ::   zdbds_mle   ! Magnitude of horizontal buoyancy gradient
381      ! Flux-gradient relationship variables
382      REAL(wp), DIMENSION(A2D(0)) ::   zshear   ! Shear production
383      !
384      REAL(wp), DIMENSION(A2D(0)) ::   zhbl_t   ! Holds boundary layer depth updated by full timestep
385      !
386      ! For calculating Ri#-dependent mixing
387      REAL(wp), DIMENSION(jpi,jpj) ::   z2du     ! u-shear^2
388      REAL(wp), DIMENSION(jpi,jpj) ::   z2dv     ! v-shear^2
389      REAL(wp)                     ::   zrimix   ! Spatial form of ri#-induced diffusion
390      !
391      ! Temporary variables
392      REAL(wp)                        ::   znd   ! Temporary non-dimensional depth
393      REAL(wp)                        ::   zz0, zz1, zfac
394      REAL(wp)                        ::   zus_x, zus_y   ! Temporary Stokes drift
395      REAL(wp), DIMENSION(A2D(0),jpk) ::   zviscos   ! Viscosity
396      REAL(wp), DIMENSION(A2D(0),jpk) ::   zdiffut   ! t-diffusivity
397      REAL(wp)                        ::   zabsstke
398      REAL(wp)                        ::   zsqrtpi, z_two_thirds, zthickness
399      REAL(wp)                        ::   z2k_times_thickness, zsqrt_depth, zexp_depth, zf, zexperfc
400      !
401      ! For debugging
402      REAL(wp), PARAMETER ::   pp_large = -1e10_wp
403      !
404      IF( ln_timing ) CALL timing_start('zdf_osm')
405      !
406      nbld(:,:)   = 0
407      nmld(:,:)   = 0
408      sustke(:,:) = pp_large
409      l_pyc(:,:)  = .FALSE.
410      l_flux(:,:) = .FALSE.
411      l_mle(:,:)  = .FALSE.
412      ! Mixed layer
413      ! No initialization of zhbl or zhml (or zdh?)
414      zhbl(:,:)     = pp_large ; zhml(:,:)     = pp_large ; zdh(:,:)      = pp_large
415      !
416      IF ( ln_osm_mle ) THEN   ! Only initialise arrays if needed
417         zdtdx(:,:)  = pp_large ; zdtdy(:,:)    = pp_large ; zdsdx(:,:)     = pp_large
418         zdsdy(:,:)  = pp_large ; dbdx_mle(:,:) = pp_large ; dbdy_mle(:,:)  = pp_large
419         zwb_fk(:,:) = pp_large ; zvel_mle(:,:) = pp_large
420         zhmle(:,:)  = pp_large ; zmld(:,:)     = pp_large
421      ENDIF
422      zhbl_t(:,:)   = pp_large
423      !
424      zdiffut(:,:,:) = 0.0_wp
425      zviscos(:,:,:) = 0.0_wp
426      !
427      ghamt(:,:,:)      = pp_large ; ghams(:,:,:)      = pp_large
428      ghamt(A2D(0),:)   = 0.0_wp   ; ghams(A2D(0),:)   = 0.0_wp
429      ghamu(:,:,:)      = pp_large ; ghamv(:,:,:)      = pp_large
430      ghamu(A2D(0),:)   = 0.0_wp   ; ghamv(A2D(0),:)   = 0.0_wp
431      !
432      zdiff_mle(:,:) = 0.0_wp
433      !
434      ! hbl = MAX(hbl,epsln)
435      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
436      ! Calculate boundary layer scales
437      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
438      !
439      ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
440      zz0 =           rn_abs   ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands
441      zz1 =  1.0_wp - rn_abs
442      DO_2D( 0, 0, 0, 0 )
443         zrad0(ji,jj)  = qsr(ji,jj) * r1_rho0_rcp   ! Surface downward irradiance (so always +ve)
444         zradh(ji,jj)  = zrad0(ji,jj) *                                &   ! Downwards irradiance at base of boundary layer
445            &            ( zz0 * EXP( -1.0_wp * hbl(ji,jj) / rn_si0 ) + zz1 * EXP( -1.0_wp * hbl(ji,jj) / rn_si1 ) )
446         zradav        = zrad0(ji,jj) *                                              &            ! Downwards irradiance averaged
447            &            ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 +   &            !    over depth of the OSBL
448            &              zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj)
449         swth0(ji,jj)  = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1)   ! Upwards surface Temperature flux for non-local term
450         swthav(ji,jj) = 0.5_wp * swth0(ji,jj) - ( 0.5_wp * ( zrad0(ji,jj) + zradh(ji,jj) ) -   &   ! Turbulent heat flux averaged
451            &                                                 zradav )                              !    over depth of OSBL
452      END_2D
453      DO_2D( 0, 0, 0, 0 )
454         sws0(ji,jj)    = -1.0_wp * ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) +   &   ! Upwards surface salinity flux
455            &                         sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)                      !    for non-local term
456         zthermal       = rab_n(ji,jj,1,jp_tem)
457         zbeta          = rab_n(ji,jj,1,jp_sal)
458         swb0(ji,jj)    = grav * zthermal * swth0(ji,jj) - grav * zbeta * sws0(ji,jj)   ! Non radiative upwards surface buoyancy flux
459         zwb0tot(ji,jj) = swb0(ji,jj) - grav * zthermal * ( zrad0(ji,jj) - zradh(ji,jj) )   ! Total upwards surface buoyancy flux
460         swsav(ji,jj)   = 0.5_wp * sws0(ji,jj)                              ! Turbulent salinity flux averaged over depth of the OBSL
461         swbav(ji,jj)   = grav  * zthermal * swthav(ji,jj) -            &   ! Turbulent buoyancy flux averaged over the depth of the
462            &             grav  * zbeta * swsav(ji,jj)                      ! OBSBL
463      END_2D
464      DO_2D( 0, 0, 0, 0 )
465         suw0(ji,jj)    = -0.5_wp * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)   ! Surface upward velocity fluxes
466         zvw0           = -0.5_wp * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
467         sustar(ji,jj)  = MAX( SQRT( SQRT( suw0(ji,jj) * suw0(ji,jj) + zvw0 * zvw0 ) ),   &   ! Friction velocity (sustar), at
468            &                  1e-8_wp )                                                      !    T-point : LMD94 eq. 2
469         scos_wind(ji,jj) = -1.0_wp * suw0(ji,jj) / ( sustar(ji,jj) * sustar(ji,jj) )
470         ssin_wind(ji,jj) = -1.0_wp * zvw0        / ( sustar(ji,jj) * sustar(ji,jj) )
471      END_2D
472      ! Calculate Stokes drift in direction of wind (sustke) and Stokes penetration depth (dstokes)
473      SELECT CASE (nn_osm_wave)
474         ! Assume constant La#=0.3
475      CASE(0)
476         DO_2D( 0, 0, 0, 0 )
477            zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2
478            zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2
479            ! Linearly
480            sustke(ji,jj)  = MAX( SQRT( zus_x * zus_x + zus_y * zus_y ), 1e-8_wp )
481            dstokes(ji,jj) = rn_osm_dstokes
482         END_2D
483         ! Assume Pierson-Moskovitz wind-wave spectrum
484      CASE(1)
485         DO_2D( 0, 0, 0, 0 )
486            ! Use wind speed wndm included in sbc_oce module
487            sustke(ji,jj)  = MAX ( 0.016_wp * wndm(ji,jj), 1e-8_wp )
488            dstokes(ji,jj) = MAX ( 0.12_wp * wndm(ji,jj)**2 / grav, 5e-1_wp )
489         END_2D
490         ! Use ECMWF wave fields as output from SBCWAVE
491      CASE(2)
492         zfac =  2.0_wp * rpi / 16.0_wp
493         !
494         DO_2D( 0, 0, 0, 0 )
495            IF ( hsw(ji,jj) > 1e-4_wp ) THEN
496               ! Use  wave fields
497               zabsstke = SQRT( ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2 )
498               sustke(ji,jj)  = MAX( ( scos_wind(ji,jj) * ut0sd(ji,jj) + ssin_wind(ji,jj)  * vt0sd(ji,jj) ), 1e-8_wp )
499               dstokes(ji,jj) = MAX( zfac * hsw(ji,jj) * hsw(ji,jj) / ( MAX( zabsstke * wmp(ji,jj), 1e-7 ) ), 5e-1_wp )
500            ELSE
501               ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run)
502               ! .. so default to Pierson-Moskowitz
503               sustke(ji,jj)  = MAX( 0.016_wp * wndm(ji,jj), 1e-8_wp )
504               dstokes(ji,jj) = MAX( 0.12_wp * wndm(ji,jj)**2 / grav, 5e-1_wp )
505            END IF
506         END_2D
507      END SELECT
508      !
509      IF (ln_zdfosm_ice_shelter) THEN
510         ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
511         DO_2D( 0, 0, 0, 0 )
512            sustke(ji,jj)  = sustke(ji,jj)  * ( 1.0_wp - fr_i(ji,jj) )
513            dstokes(ji,jj) = dstokes(ji,jj) * ( 1.0_wp - fr_i(ji,jj) )
514         END_2D
515      END IF
516      !
517      SELECT CASE (nn_osm_SD_reduce)
518         ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van Roekel (2012) or Grant (2020).
519      CASE(0)
520         ! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas.
521         ! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation.
522         ! It could represent the effects of the spread of wave directions around the mean wind. The effect of this adjustment needs to be tested.
523         IF(nn_osm_wave > 0) THEN
524            sustke(A2D(0)) = rn_zdfosm_adjust_sd * sustke(A2D(0))
525         END IF
526      CASE(1)
527         ! Van Roekel (2012): consider average SD over top 10% of boundary layer
528         ! Assumes approximate depth profile of SD from Breivik (2016)
529         zsqrtpi = SQRT(rpi)
530         z_two_thirds = 2.0_wp / 3.0_wp
531         DO_2D( 0, 0, 0, 0 )
532            zthickness = rn_osm_hblfrac*hbl(ji,jj)
533            z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp )
534            zsqrt_depth = SQRT( z2k_times_thickness )
535            zexp_depth  = EXP( -1.0_wp * z2k_times_thickness )
536            sustke(ji,jj) = sustke(ji,jj) * ( 1.0_wp - zexp_depth -   &
537               &                              z_two_thirds * ( zsqrtpi * zsqrt_depth * z2k_times_thickness * ERFC(zsqrt_depth) +   &
538               &                                               1.0_wp - ( 1.0_wp + z2k_times_thickness ) * zexp_depth ) ) /        &
539               &            z2k_times_thickness
540         END_2D
541      CASE(2)
542         ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer
543         ! Assumes approximate depth profile of SD from Breivik (2016)
544         zsqrtpi = SQRT(rpi)
545         DO_2D( 0, 0, 0, 0 )
546            zthickness = rn_osm_hblfrac*hbl(ji,jj)
547            z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp )
548            IF( z2k_times_thickness < 50.0_wp ) THEN
549               zsqrt_depth = SQRT( z2k_times_thickness )
550               zexperfc    = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP( z2k_times_thickness )
551            ELSE
552               ! Asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness
553               ! See Abramowitz and Stegun, Eq. 7.1.23
554               ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3)
555               zexperfc = ( ( -1.875_wp / z2k_times_thickness + 0.75_wp ) / z2k_times_thickness - 0.5_wp ) /   &
556                  &       z2k_times_thickness + 1.0_wp
557            END IF
558            zf = z2k_times_thickness * ( 1.0_wp / zexperfc - 1.0_wp )
559            dstokes(ji,jj) = 5.97_wp * zf * dstokes(ji,jj)
560            sustke(ji,jj)  = sustke(ji,jj) * EXP( z2k_times_thickness * ( 1.0_wp / ( 2.0_wp * zf ) - 1.0_wp ) ) *   &
561               &             ( 1.0_wp - zexperfc )
562         END_2D
563      END SELECT
564      !
565      ! Langmuir velocity scale (swstrl), La # (sla)
566      ! Mixed scale (svstr), convective velocity scale (swstrc)
567      DO_2D( 0, 0, 0, 0 )
568         ! Langmuir velocity scale (swstrl), at T-point
569         swstrl(ji,jj) = ( sustar(ji,jj) * sustar(ji,jj) * sustke(ji,jj) )**pthird
570         sla(ji,jj)    = MAX( MIN( SQRT( sustar(ji,jj) / ( swstrl(ji,jj) + epsln ) )**3, 4.0_wp ), 0.2_wp )
571         IF ( sla(ji,jj) > 0.45_wp ) dstokes(ji,jj) = MIN( dstokes(ji,jj), 0.5_wp * hbl(ji,jj) )
572         ! Velocity scale that tends to sustar for large Langmuir numbers
573         svstr(ji,jj) = ( swstrl(ji,jj)**3 + ( 1.0_wp - EXP( -0.5_wp * sla(ji,jj)**2 ) ) * sustar(ji,jj) * sustar(ji,jj) *   &
574            &                                sustar(ji,jj) )**pthird
575         !
576         ! Limit maximum value of Langmuir number as approximate treatment for shear turbulence
577         ! Note sustke and swstrl are not amended
578         !
579         ! Get convective velocity (swstrc), stabilty scale (shol) and logical conection flag l_conv
580         IF ( swbav(ji,jj) > 0.0_wp ) THEN
581            swstrc(ji,jj) = ( 2.0_wp * swbav(ji,jj) * 0.9_wp * hbl(ji,jj) )**pthird
582            shol(ji,jj)   = -0.9_wp * hbl(ji,jj) * 2.0_wp * swbav(ji,jj) / ( svstr(ji,jj)**3 + epsln )
583         ELSE
584            swstrc(ji,jj) = 0.0_wp
585            shol(ji,jj)   = -1.0_wp * hbl(ji,jj) * 2.0_wp * swbav(ji,jj) / ( svstr(ji,jj)**3  + epsln )
586         ENDIF
587      END_2D
588      !
589      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
590      ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
591      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
592      ! BL must be always 4 levels deep.
593      ! For calculation of lateral buoyancy gradients for FK in
594      ! zdf_osm_zmld_horizontal_gradients need halo values for nbld, so must
595      ! previously exist for hbl also.
596      !
597      ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway
598      ! ##########################################################################
599      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) )
600      nbld(:,:) = 4
601      DO_3D( 1, 1, 1, 1, 5, jpkm1 )
602         IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
603            nbld(ji,jj) = MIN(mbkt(ji,jj)-2, jk)
604         ENDIF
605      END_3D
606      ! ##########################################################################
607      !
608      DO_2D( 0, 0, 0, 0 )
609         zhbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm)
610         nmld(ji,jj) = MAX( 3, nbld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji,jj,nbld(ji,jj)-1,Kmm) ), 1 ) )
611         zhml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm)
612         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
613      END_2D
614      !
615      ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere
616      jp_ext(:,:) = 1   ! ag 19/03
617      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      &
618         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   &
619         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
620      jp_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03
621      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          &
622         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   &
623         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )
624      ! Velocity components in frame aligned with surface stress
625      CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml  )
626      CALL zdf_osm_velocity_rotation( av_du_ml, av_dv_ml )
627      CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl  )
628      CALL zdf_osm_velocity_rotation( av_du_bl, av_dv_bl )
629      !
630      ! Determine the state of the OSBL, stable/unstable, shear/no shear
631      CALL zdf_osm_osbl_state( Kmm, zwb_ent, zwb_min, zshear, zhbl,     &
632         &                     zhml, zdh )
633      !
634      IF ( ln_osm_mle ) THEN
635         ! Fox-Kemper Scheme
636         mld_prof = 4
637         DO_3D( 0, 0, 0, 0, 5, jpkm1 )
638            IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
639         END_3D
640         CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D(0)), av_t_mle, av_s_mle,   &
641            &                           av_b_mle, av_u_mle, av_v_mle )
642         !
643         DO_2D( 0, 0, 0, 0 )
644            zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
645         END_2D
646         !
647         ! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients
648         CALL zdf_osm_zmld_horizontal_gradients( Kmm, zmld, zdtdx, zdtdy, zdsdx,   &
649            &                                    zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
650         ! Calculate max vertical FK flux zwb_fk & set logical descriptors
651         CALL zdf_osm_osbl_state_fk( Kmm, zwb_fk, zhbl, zhmle, zwb_ent,   &
652            &                        zdbds_mle )
653         ! Recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle
654         CALL zdf_osm_mle_parameters( Kmm, mld_prof, zmld, zhmle, zvel_mle,   &
655            &                         zdiff_mle, zdbds_mle, zhbl, zwb0tot )
656      ELSE    ! ln_osm_mle
657         ! FK not selected, Boundary Layer only.
658         l_pyc(:,:)  = .TRUE.
659         l_flux(:,:) = .FALSE.
660         l_mle(:,:)  = .FALSE.
661         DO_2D( 0, 0, 0, 0 )
662            IF ( l_conv(ji,jj) .AND. av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE.
663         END_2D
664      ENDIF   ! ln_osm_mle
665      !
666      !! External gradient below BL needed both with and w/o FK
667      CALL zdf_osm_external_gradients( Kmm, nbld(A2D(0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03
668      !
669      ! Test if pycnocline well resolved
670      !      DO_2D( 0, 0, 0, 0 )                                         Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity
671      !         IF (l_conv(ji,jj) ) THEN                                  should account for this.
672      !            ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,nbld(ji,jj),Kmm)
673      !            IF ( ztmp > 6 ) THEN
674      !               ! pycnocline well resolved
675      !               jp_ext(ji,jj) = 1
676      !            ELSE
677      !               ! pycnocline poorly resolved
678      !               jp_ext(ji,jj) = 0
679      !            ENDIF
680      !         ELSE
681      !            ! Stable conditions
682      !            jp_ext(ji,jj) = 0
683      !         ENDIF
684      !      END_2D
685      !
686      ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here..
687      jp_ext(:,:) = 1   ! ag 19/03
688      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      &
689         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   &
690         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
691      jp_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03
692      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          &
693         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   &
694         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )   ! ag 19/03
695      !
696      ! Rate of change of hbl
697      CALL zdf_osm_calculate_dhdt( zdhdt, zhbl, zdh, zwb_ent, zwb_min,   &
698         &                         zdbdz_bl_ext, zwb_fk_b, zwb_fk, zvel_mle )
699      ! Test if surface boundary layer coupled to bottom
700      l_coup(:,:) = .FALSE.   ! ag 19/03
701      DO_2D( 0, 0, 0, 0 )
702         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
703         ! Adjustment to represent limiting by ocean bottom
704         IF ( mbkt(ji,jj) > 2 ) THEN   ! To ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access
705            IF ( zhbl_t(ji,jj) > gdepw(ji, jj,mbkt(ji,jj)-2,Kmm) ) THEN
706               zhbl_t(ji,jj) = MIN( zhbl_t(ji,jj), gdepw(ji,jj,mbkt(ji,jj)-2,Kmm) )   ! ht(:,:))
707               l_pyc(ji,jj)  = .FALSE.
708               l_coup(ji,jj) = .TRUE.   ! ag 19/03
709            END IF
710         END IF
711      END_2D
712      !
713      nmld(:,:) = nbld(A2D(0))           ! use nmld to hold previous blayer index
714      nbld(:,:) = 4
715      !
716      DO_3D( 0, 0, 0, 0, 4, jpkm1 )
717         IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
718            nbld(ji,jj) = jk
719         END IF
720      END_3D
721      !
722      !
723      ! Step through model levels taking account of buoyancy change to determine the effect on dhdt
724      !
725      CALL zdf_osm_timestep_hbl( Kmm, zdhdt, zhbl, zhbl_t, zwb_ent,   &
726         &                       zwb_fk_b )
727      ! Is external level in bounds?
728      !
729      ! Recalculate BL averages and differences using new BL depth
730      jp_ext(:,:) = 1   ! ag 19/03
731      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      &
732         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   &
733         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
734      !
735      CALL zdf_osm_pycnocline_thickness( Kmm, zdh, zhml, zdhdt, zhbl,   &
736         &                               zwb_ent, zdbdz_bl_ext, zwb_fk_b )
737      !
738      ! Reset l_pyc before calculating terms in the flux-gradient relationship
739      DO_2D( 0, 0, 0, 0 )
740         IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR.   &
741            & nbld(ji,jj) - nmld(ji,jj) == 1   .OR. zdhdt(ji,jj) < 0.0_wp ) THEN   ! ag 19/03
742            l_pyc(ji,jj) = .FALSE.   ! ag 19/03
743            IF ( nbld(ji,jj) >= mbkt(ji,jj) -2 ) THEN
744               nmld(ji,jj) = nbld(ji,jj) - 1                                               ! ag 19/03
745               zdh(ji,jj)  = gdepw(ji,jj,nbld(ji,jj),Kmm) - gdepw(ji,jj,nmld(ji,jj),Kmm)   ! ag 19/03
746               zhml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm)                                  ! ag 19/03
747               dh(ji,jj)   = zdh(ji,jj)                                                    ! ag 19/03 
748               hml(ji,jj)  = hbl(ji,jj) - dh(ji,jj)                                        ! ag 19/03
749            ENDIF
750         ENDIF                                              ! ag 19/03
751      END_2D
752      !
753      dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D(0))/ 3.0_wp )   ! Limit delta for shallow boundary layers for calculating
754      !                                                       !    flux-gradient terms
755      !
756      ! Average over the depth of the mixed layer in the convective boundary layer
757      !      jp_ext = nbld - nmld + 1
758      ! Recalculate ML averages and differences using new ML depth
759      jp_ext(:,:) = nbld(A2D(0)) - nmld(A2D(0)) + jp_ext(:,:) + 1   ! ag 19/03
760      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,    &
761         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   &
762         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )
763      !
764      CALL zdf_osm_external_gradients( Kmm, nbld(A2D(0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
765      ! Rotate mean currents and changes onto wind aligned co-ordinates
766      CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml )
767      CALL zdf_osm_velocity_rotation( av_du_ml, av_dv_ml )
768      CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl )
769      CALL zdf_osm_velocity_rotation( av_du_bl, av_dv_bl )
770      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
771      ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
772      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
773      CALL zdf_osm_diffusivity_viscosity( Kbb, Kmm, zdiffut, zviscos, zhbl,    &
774         &                                zhml, zdh, zdhdt, zshear, zwb_ent,   &
775         &                                zwb_min )
776      !
777      ! Calculate non-gradient components of the flux-gradient relationships
778      ! --------------------------------------------------------------------
779      jp_ext(:,:) = 1   ! ag 19/03
780      CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh,                              &
781         &                    zdhdt, zshear, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext,   &
782         &                    zdiffut, zviscos )
783      !
784      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
785      ! Need to put in code for contributions that are applied explicitly to
786      ! the prognostic variables
787      !  1. Entrainment flux
788      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
789      !
790      ! Rotate non-gradient velocity terms back to model reference frame
791      CALL zdf_osm_velocity_rotation( ghamu(A2D(0),:), ghamv(A2D(0),:), .FALSE., 2, nbld(A2D(0)) )
792      !
793      ! KPP-style Ri# mixing
794      IF ( ln_kpprimix ) THEN
795         jkflt = jpk
796         DO_2D( 0, 0, 0, 0 )
797            IF ( nbld(ji,jj) < jkflt ) jkflt = nbld(ji,jj)
798         END_2D
799         DO jk = jkflt+1, jpkm1
800            ! Shear production at uw- and vw-points (energy conserving form)
801            DO_2D( 1, 0, 1, 0 )
802               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) ) *   &
803                  &          wumask(ji,jj,jk) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) )
804               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) ) *   &
805                  &          wvmask(ji,jj,jk) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) )
806            END_2D
807            DO_2D( 0, 0, 0, 0 )
808               IF ( jk > nbld(ji,jj) ) THEN
809                  ! Shear prod. at w-point weightened by mask
810                  zesh2 = ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) +   &
811                     &    ( z2dv(ji,jj-1) + z2dv(ji,jj) ) / MAX( 1.0_wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
812                  ! Local Richardson number
813                  zri     = MAX( rn2b(ji,jj,jk), 0.0_wp ) / MAX( zesh2, epsln )
814                  zfri    = MIN( zri / rn_riinfty, 1.0_wp )
815                  zfri    = ( 1.0_wp - zfri * zfri )
816                  zrimix  =  zfri * zfri  * zfri * wmask(ji, jj, jk)
817                  zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zrimix*rn_difri )
818                  zviscos(ji,jj,jk) = MAX( zviscos(ji,jj,jk), zrimix*rn_difri )
819               END IF
820            END_2D
821         END DO
822      END IF   ! ln_kpprimix = .true.
823      !
824      ! KPP-style set diffusivity large if unstable below BL
825      IF ( ln_convmix) THEN
826         DO_2D( 0, 0, 0, 0 )
827            DO jk = nbld(ji,jj) + 1, jpkm1
828               IF ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1e-12_wp ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) )
829            END DO
830         END_2D
831      END IF   ! ln_convmix = .true.
832      !
833      IF ( ln_osm_mle ) THEN   ! Set up diffusivity and non-gradient mixing
834         DO_2D( 0, 0, 0, 0 )
835            IF ( l_flux(ji,jj) ) THEN   ! MLE mixing extends below boundary layer
836               ! Calculate MLE flux contribution from surface fluxes
837               DO jk = 1, nbld(ji,jj)
838                  znd = gdepw(ji,jj,jk,Kmm) / MAX( zhbl(ji,jj), epsln )
839                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0_wp - znd )
840                  ghams(ji,jj,jk) = ghams(ji,jj,jk) - sws0(ji,jj) * ( 1.0_wp - znd )
841               END DO
842               DO jk = 1, mld_prof(ji,jj)
843                  znd = gdepw(ji,jj,jk,Kmm) / MAX( zhmle(ji,jj), epsln )
844                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( swth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0_wp - znd )
845                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + sws0(ji,jj) * ( 1.0_wp -znd )
846               END DO
847               ! Viscosity for MLEs
848               DO jk = 1, mld_prof(ji,jj)
849                  znd = -1.0_wp * gdepw(ji,jj,jk,Kmm) / MAX( zhmle(ji,jj), epsln )
850                  zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0_wp - ( 2.0_wp * znd + 1.0_wp )**2 ) *   &
851                     &                                    ( 1.0_wp + 5.0_wp / 21.0_wp * ( 2.0_wp * znd + 1.0_wp )**2 )
852               END DO
853            ELSE   ! Surface transports limited to OSBL
854               ! Viscosity for MLEs
855               DO jk = 1, mld_prof(ji,jj)
856                  znd = -1.0_wp * gdepw(ji,jj,jk,Kmm) / MAX( zhmle(ji,jj), epsln )
857                  zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0_wp - ( 2.0_wp * znd + 1.0_wp )**2 ) *   &
858                     &                                    ( 1.0_wp + 5.0_wp / 21.0_wp * ( 2.0_wp * znd + 1.0_wp )**2 )
859               END DO
860            END IF
861         END_2D
862      ENDIF
863      !
864      ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
865      ! CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp )
866      ! GN 25/8: need to change tmask --> wmask
867      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
868         p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
869         p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
870      END_3D
871      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and
872      !    v grids
873      CALL lbc_lnk( 'zdfosm', p_avt, 'W', 1.0_wp,   &
874         &                    p_avm, 'W', 1.0_wp,   &
875         &                    ghamu, 'W', 1.0_wp,   &
876         &                    ghamv, 'W', 1.0_wp )
877      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
878         ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji+1,jj,jk) ) *   &
879            &              umask(ji,jj,jk)
880         ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) *   &
881            &              vmask(ji,jj,jk)
882         ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk)
883         ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk)
884      END_3D
885      ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)
886      CALL lbc_lnk( 'zdfosm', hbl,  'T', 1.0_wp,   &
887         &                    dh,   'T', 1.0_wp,   &
888         &                    hmle, 'T', 1.0_wp )
889      ! Lateral boundary conditions on final outputs for gham[ts], on W-grid  (sign unchanged)
890      ! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid  (sign changed)
891      CALL lbc_lnk( 'zdfosm', ghamt, 'W',  1.0_wp,   &
892         &                    ghams, 'W',  1.0_wp,   &
893         &                    ghamu, 'U', -1.0_wp,   &
894         &                    ghamv, 'V', -1.0_wp )
895      !
896      IF ( ln_dia_osm ) THEN
897         SELECT CASE (nn_osm_wave)
898            ! Stokes drift set by assumimg onstant La#=0.3 (=0) or Pierson-Moskovitz spectrum (=1)
899         CASE(0:1)
900            IF ( iom_use("us_x") ) THEN                                                           ! x surface Stokes drift
901               osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke * scos_wind
902               CALL iom_put( "us_x", osmdia2d )
903            END IF
904            IF ( iom_use("us_y") ) THEN                                                           ! y surface Stokes drift
905               osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke * ssin_wind
906               CALL iom_put( "us_y", osmdia2d )
907            END IF
908            IF ( iom_use("wind_wave_abs_power") ) THEN
909               osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke
910               CALL iom_put( "wind_wave_abs_power", osmdia2d )
911            END IF
912            ! Stokes drift read in from sbcwave  (=2).
913         CASE(2:3)
914            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )                     ! x surface Stokes drift
915            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )                     ! y surface Stokes drift
916            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                         ! Wave mean period
917            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                         ! Significant wave height
918            IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", ( 2.0_wp * rpi * 1.026_wp /      &   ! Wave mean period from NP
919               &                                               ( 0.877_wp * grav ) ) *        &   !    spectrum
920               &                                               wndm * tmask(:,:,1) )
921            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", ( 0.22_wp / grav ) * wndm**2 *   &   ! Significant wave
922               &                                             tmask(:,:,1) )                       !    height from NP spectrum
923            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                      ! U_10
924            IF ( iom_use("wind_wave_abs_power") ) THEN
925               osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(:,:,1) * sustar**2 * SQRT( ut0sd**2 + vt0sd**2 )
926               CALL iom_put( "wind_wave_abs_power", osmdia2d )
927            END IF
928         END SELECT
929         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )                      ! <Tw_NL>
930         IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams )                      ! <Sw_NL>
931         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )                      ! <uw_NL>
932         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )                      ! <vw_NL>
933         IF ( iom_use("zwth0") ) THEN                                                      ! <Tw_0>
934            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swth0;     CALL iom_put( "zwth0",     osmdia2d )
935         END IF
936         IF ( iom_use("zws0") ) THEN                                                       ! <Sw_0>
937            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sws0;      CALL iom_put( "zws0",      osmdia2d )
938         END IF
939         IF ( iom_use("zwb0") ) THEN                                                       ! <Sw_0>
940            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swb0;      CALL iom_put( "zwb0",      osmdia2d )
941         END IF
942         IF ( iom_use("zwbav") ) THEN                                                      ! Upward BL-avged turb buoyancy flux
943            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swth0;     CALL iom_put( "zwbav",     osmdia2d )
944         END IF
945         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                     ! Boundary-layer depth
946         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld )                  ! Boundary-layer max k
947         IF ( iom_use("zdt_bl") ) THEN                                                     ! dt at ml base
948            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dt_bl;  CALL iom_put( "zdt_bl", osmdia2d )
949         END IF
950         IF ( iom_use("zds_bl") ) THEN                                                     ! ds at ml base
951            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_ds_bl;  CALL iom_put( "zds_bl", osmdia2d )
952         END IF
953         IF ( iom_use("zdb_bl") ) THEN                                                     ! db at ml base
954            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_db_bl;  CALL iom_put( "zdb_bl", osmdia2d )
955         END IF
956         IF ( iom_use("zdu_bl") ) THEN                                                     ! du at ml base
957            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_du_bl;  CALL iom_put( "zdu_bl", osmdia2d )
958         END IF
959         IF ( iom_use("zdv_bl") ) THEN                                                     ! dv at ml base
960            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dv_bl;  CALL iom_put( "zdv_bl", osmdia2d )
961         END IF
962         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )                        ! Initial boundary-layer depth
963         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )                     ! Initial boundary-layer depth
964         IF ( iom_use("zdt_ml") ) THEN                                                     ! dt at ml base
965            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dt_ml;  CALL iom_put( "zdt_ml", osmdia2d )
966         END IF
967         IF ( iom_use("zds_ml") ) THEN                                                     ! ds at ml base
968            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_ds_ml;  CALL iom_put( "zds_ml", osmdia2d )
969         END IF
970         IF ( iom_use("zdb_ml") ) THEN                                                     ! db at ml base
971            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_db_ml;  CALL iom_put( "zdb_ml", osmdia2d )
972         END IF
973         IF ( iom_use("dstokes") ) THEN                                                    ! Stokes drift penetration depth
974            osmdia2d(A2D(0)) = tmask(A2D(0),1) * dstokes;   CALL iom_put( "dstokes",   osmdia2d )
975         END IF
976         IF ( iom_use("zustke") ) THEN                                                     ! Stokes drift magnitude at T-points
977            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke;    CALL iom_put( "zustke",    osmdia2d )
978         END IF
979         IF ( iom_use("zwstrc") ) THEN                                                     ! Convective velocity scale
980            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swstrc;    CALL iom_put( "zwstrc",    osmdia2d )
981         END IF
982         IF ( iom_use("zwstrl") ) THEN                                                     ! Langmuir velocity scale
983            osmdia2d(A2D(0)) = tmask(A2D(0),1) * swstrl;    CALL iom_put( "zwstrl",    osmdia2d )
984         END IF
985         IF ( iom_use("zustar") ) THEN                                                     ! Friction velocity scale
986            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustar;    CALL iom_put( "zustar",    osmdia2d )
987         END IF
988         IF ( iom_use("zvstr") ) THEN                                                      ! Mixed velocity scale
989            osmdia2d(A2D(0)) = tmask(A2D(0),1) * svstr;     CALL iom_put( "zvstr",     osmdia2d )
990         END IF
991         IF ( iom_use("zla") ) THEN                                                        ! Langmuir #
992            osmdia2d(A2D(0)) = tmask(A2D(0),1) * sla;       CALL iom_put( "zla",       osmdia2d )
993         END IF
994         IF ( iom_use("wind_power") ) THEN                                                 ! BL depth internal to zdf_osm routine
995            osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**3
996            CALL iom_put( "wind_power", osmdia2d )
997         END IF
998         IF ( iom_use("wind_wave_power") ) THEN
999            osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke
1000            CALL iom_put( "wind_wave_power", osmdia2d )
1001         END IF
1002         IF ( iom_use("zhbl") ) THEN                                                       ! BL depth internal to zdf_osm routine
1003            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zhbl;      CALL iom_put( "zhbl",      osmdia2d )
1004         END IF
1005         IF ( iom_use("zhml") ) THEN                                                       ! ML depth internal to zdf_osm routine
1006            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zhml;      CALL iom_put( "zhml",      osmdia2d )
1007         END IF
1008         IF ( iom_use("imld") ) THEN                                                       ! Index for ML depth internal to zdf_osm routine
1009            osmdia2d(A2D(0)) = tmask(A2D(0),1) * nmld;      CALL iom_put( "imld",      osmdia2d )
1010         END IF
1011         IF ( iom_use("jp_ext") ) THEN                                                     ! =1 if pycnocline resolved internal to zdf_osm routine
1012            osmdia2d(A2D(0)) = tmask(A2D(0),1) * jp_ext;    CALL iom_put( "jp_ext",    osmdia2d )
1013         END IF
1014         IF ( iom_use("j_ddh") ) THEN                                                      ! Index forpyc thicknessh internal to zdf_osm routine
1015            osmdia2d(A2D(0)) = tmask(A2D(0),1) * n_ddh;     CALL iom_put( "j_ddh",     osmdia2d )
1016         END IF
1017         IF ( iom_use("zshear") ) THEN                                                     ! Shear production of TKE internal to zdf_osm routine
1018            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zshear;    CALL iom_put( "zshear",    osmdia2d )
1019         END IF
1020         IF ( iom_use("zdh") ) THEN                                                        ! Pyc thicknessh internal to zdf_osm routine
1021            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdh;       CALL iom_put( "zdh",       osmdia2d )
1022         END IF
1023         IF ( iom_use("zhol") ) THEN                                                       ! ML depth internal to zdf_osm routine
1024            osmdia2d(A2D(0)) = tmask(A2D(0),1) * shol;      CALL iom_put( "zhol",      osmdia2d )
1025         END IF
1026         IF ( iom_use("zwb_ent") ) THEN                                                    ! Upward turb buoyancy entrainment flux
1027            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_ent;   CALL iom_put( "zwb_ent",   osmdia2d )
1028         END IF
1029         IF ( iom_use("zt_ml") ) THEN                                                      ! Average T in ML
1030            osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_t_ml;   CALL iom_put( "zt_ml",     osmdia2d )
1031         END IF
1032         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )                  ! FK layer depth
1033         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )                  ! FK target layer depth
1034         IF ( iom_use("zwb_fk") ) THEN                                                     ! FK b flux
1035            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_fk;    CALL iom_put( "zwb_fk",    osmdia2d )
1036         END IF
1037         IF ( iom_use("zwb_fk_b") ) THEN                                                   ! FK b flux averaged over ML
1038            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_fk_b;  CALL iom_put( "zwb_fk_b",  osmdia2d )
1039         END IF
1040         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )      ! FK layer max k
1041         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )               ! FK dtdx at u-pt
1042         IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy )               ! FK dtdy at v-pt
1043         IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx )               ! FK dtdx at u-pt
1044         IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy )               ! FK dsdy at v-pt
1045         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )      ! FK dbdx at u-pt
1046         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )      ! FK dbdy at v-pt
1047         IF ( iom_use("zdiff_mle") ) THEN                                                  ! FK diff in MLE at t-pt
1048            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d )
1049         END IF
1050         IF ( iom_use("zvel_mle") ) THEN                                                   ! FK diff in MLE at t-pt
1051            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zvel_mle",  osmdia2d )
1052         END IF
1053      END IF
1054      IF( ln_timing ) CALL timing_stop('zdf_osm')
1055
1056   END SUBROUTINE zdf_osm
1057
1058   SUBROUTINE zdf_osm_vertical_average( Kbb, Kmm, knlev, pt, ps,   &
1059      &                                 pb, pu, pv, kp_ext, pdt,   &
1060      &                                 pds, pdb, pdu, pdv )
1061      !!---------------------------------------------------------------------
1062      !!                ***  ROUTINE zdf_vertical_average  ***
1063      !!
1064      !! ** Purpose : Determines vertical averages from surface to knlev,
1065      !!              and optionally the differences between these vertical
1066      !!              averages and values at an external level
1067      !!
1068      !! ** Method  : Averages are calculated from the surface to knlev.
1069      !!              The external level used to calculate differences is
1070      !!              knlev+kp_ext
1071      !!----------------------------------------------------------------------
1072      INTEGER,                     INTENT(in   )           ::   Kbb, Kmm   ! Ocean time-level indices
1073      INTEGER,  DIMENSION(A2D(0)), INTENT(in   )           ::   knlev      ! Number of levels to average over.
1074      REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pt, ps     ! Average temperature and salinity
1075      REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pb         ! Average buoyancy
1076      REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pu, pv     ! Average current components
1077      INTEGER,  DIMENSION(A2D(0)), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets
1078      REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdt, pds   ! Difference between average temperature, salinity,
1079      REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdb        ! buoyancy,
1080      REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdu, pdv   ! velocity components and the OSBL
1081      !
1082      INTEGER                     ::   jk, jkflt, jkmax, ji, jj   ! Loop indices
1083      INTEGER                     ::   ibld_ext                   ! External-layer index
1084      REAL(wp), DIMENSION(A2D(0)) ::   zthick                     ! Layer thickness
1085      REAL(wp)                    ::   zthermal, zbeta            ! Thermal/haline expansion/contraction coefficients
1086      !!----------------------------------------------------------------------
1087      !
1088      IF( ln_timing ) CALL timing_start('zdf_osm_va')
1089      !
1090      ! Averages over depth of boundary layer
1091      pt(:,:)     = 0.0_wp
1092      ps(:,:)     = 0.0_wp
1093      pu(:,:)     = 0.0_wp
1094      pv(:,:)     = 0.0_wp
1095      zthick(:,:) = epsln
1096      jkflt = jpk
1097      jkmax = 0
1098      DO_2D( 0, 0, 0, 0 )
1099         IF ( knlev(ji,jj) < jkflt ) jkflt = knlev(ji,jj)
1100         IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj)
1101      END_2D
1102      DO_3D( 0, 0, 0, 0, 2, jkflt )   ! Upper, flat part of layer
1103         zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
1104         pt(ji,jj)     = pt(ji,jj)     + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
1105         ps(ji,jj)     = ps(ji,jj)     + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
1106         pu(ji,jj)     = pu(ji,jj)     + e3t(ji,jj,jk,Kmm) *                                        &
1107            &                               ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) /           &
1108            &                               MAX( 1.0_wp , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
1109         pv(ji,jj)     = pv(ji,jj)     + e3t(ji,jj,jk,Kmm) *                                        &
1110            &                               ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) /           &
1111            &                               MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )         
1112      END_3D
1113      DO_3D( 0, 0, 0, 0, jkflt+1, jkmax )   ! Lower, non-flat part of layer
1114         IF ( knlev(ji,jj) >= jk ) THEN
1115            zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
1116            pt(ji,jj)     = pt(ji,jj)     + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
1117            ps(ji,jj)     = ps(ji,jj)     + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
1118            pu(ji,jj)     = pu(ji,jj)     + e3t(ji,jj,jk,Kmm) *                                        &
1119               &                               ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) /           &
1120               &                               MAX( 1.0_wp , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
1121            pv(ji,jj)     = pv(ji,jj)     + e3t(ji,jj,jk,Kmm) *                                        &
1122               &                               ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) /           &
1123               &                               MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
1124         END IF
1125      END_3D
1126      DO_2D( 0, 0, 0, 0 )
1127         pt(ji,jj) = pt(ji,jj) / zthick(ji,jj)
1128         ps(ji,jj) = ps(ji,jj) / zthick(ji,jj)
1129         pu(ji,jj) = pu(ji,jj) / zthick(ji,jj)
1130         pv(ji,jj) = pv(ji,jj) / zthick(ji,jj)
1131         zthermal  = rab_n(ji,jj,1,jp_tem)   ! ideally use nbld not 1??
1132         zbeta     = rab_n(ji,jj,1,jp_sal)
1133         pb(ji,jj) = grav * zthermal * pt(ji,jj) - grav * zbeta * ps(ji,jj)
1134      END_2D
1135      !
1136      ! Differences between vertical averages and values at an external layer
1137      IF ( PRESENT( kp_ext ) ) THEN
1138         DO_2D( 0, 0, 0, 0 )
1139            ibld_ext = knlev(ji,jj) + kp_ext(ji,jj)
1140            IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN   ! ag 09/03
1141               ! Two external levels are available
1142               pdt(ji,jj) = pt(ji,jj) - ts(ji,jj,ibld_ext,jp_tem,Kmm)
1143               pds(ji,jj) = ps(ji,jj) - ts(ji,jj,ibld_ext,jp_sal,Kmm)
1144               pdu(ji,jj) = pu(ji,jj) - ( uu(ji,jj,ibld_ext,Kbb) + uu(ji-1,jj,ibld_ext,Kbb ) ) /              &
1145                  &                        MAX(1.0_wp , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) )
1146               pdv(ji,jj) = pv(ji,jj) - ( vv(ji,jj,ibld_ext,Kbb) + vv(ji,jj-1,ibld_ext,Kbb ) ) /              &
1147                  &                        MAX(1.0_wp , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) )
1148               zthermal   = rab_n(ji,jj,1,jp_tem)   ! ideally use nbld not 1??
1149               zbeta      = rab_n(ji,jj,1,jp_sal)
1150               pdb(ji,jj) = grav * zthermal * pdt(ji,jj) - grav * zbeta * pds(ji,jj)
1151            ELSE
1152               pdt(ji,jj) = 0.0_wp
1153               pds(ji,jj) = 0.0_wp
1154               pdu(ji,jj) = 0.0_wp
1155               pdv(ji,jj) = 0.0_wp
1156               pdb(ji,jj) = 0.0_wp
1157            ENDIF
1158         END_2D
1159      END IF
1160      !
1161      IF( ln_timing ) CALL timing_stop('zdf_osm_va')
1162      !
1163   END SUBROUTINE zdf_osm_vertical_average
1164
1165   SUBROUTINE zdf_osm_velocity_rotation_2d( pu, pv, fwd )
1166      !!---------------------------------------------------------------------
1167      !!            ***  ROUTINE zdf_velocity_rotation_2d  ***
1168      !!
1169      !! ** Purpose : Rotates frame of reference of velocity components pu and
1170      !!              pv (2d)
1171      !!
1172      !! ** Method : The velocity components are rotated into (fwd=.TRUE.) or
1173      !!             from (fwd=.FALSE.) the frame specified by scos_wind and
1174      !!             ssin_wind
1175      !!
1176      !!----------------------------------------------------------------------     
1177      REAL(wp),           INTENT(inout), DIMENSION(A2D(0)) ::   pu, pv           ! Components of current
1178      LOGICAL,  OPTIONAL, INTENT(in   )                    ::   fwd              ! Forward (default) or reverse rotation
1179      !
1180      INTEGER  ::   ji, jj       ! Loop indices
1181      REAL(wp) ::   ztmp, zfwd   ! Auxiliary variables
1182      !
1183      IF( ln_timing ) CALL timing_start('zdf_osm_vr')
1184      !
1185      zfwd = 1.0_wp
1186      IF( PRESENT(fwd) .AND. ( fwd==.FALSE. ) ) zfwd = -1.0_wp
1187      DO_2D( 0, 0, 0, 0 )
1188         ztmp      = pu(ji,jj)
1189         pu(ji,jj) = pu(ji,jj) * scos_wind(ji,jj) + zfwd * pv(ji,jj) * ssin_wind(ji,jj)
1190         pv(ji,jj) = pv(ji,jj) * scos_wind(ji,jj) - zfwd * ztmp      * ssin_wind(ji,jj)
1191      END_2D
1192      !
1193      IF( ln_timing ) CALL timing_stop('zdf_osm_vr')
1194      !
1195   END SUBROUTINE zdf_osm_velocity_rotation_2d
1196
1197   SUBROUTINE zdf_osm_velocity_rotation_3d( pu, pv, fwd, ktop, knlev )
1198      !!---------------------------------------------------------------------
1199      !!            ***  ROUTINE zdf_velocity_rotation_3d  ***
1200      !!
1201      !! ** Purpose : Rotates frame of reference of velocity components pu and
1202      !!              pv (3d)
1203      !!
1204      !! ** Method : The velocity components are rotated into (fwd=.TRUE.) or
1205      !!             from (fwd=.FALSE.) the frame specified by scos_wind and
1206      !!             ssin_wind; optionally, the rotation can be restricted at
1207      !!             each water column to span from the a minimum index ktop to
1208      !!             the depth index specified in array knlev
1209      !!
1210      !!----------------------------------------------------------------------     
1211      REAL(wp),           INTENT(inout), DIMENSION(A2D(0),jpk) ::   pu, pv           ! Components of current
1212      LOGICAL,  OPTIONAL, INTENT(in   )                        ::   fwd              ! Forward (default) or reverse rotation
1213      INTEGER,  OPTIONAL, INTENT(in   )                        ::   ktop             ! Minimum depth index
1214      INTEGER,  OPTIONAL, INTENT(in   ), DIMENSION(A2D(0))     ::   knlev            ! Array of maximum depth indices
1215      !
1216      INTEGER  ::   ji, jj, jk, jktop, jkmax   ! Loop indices
1217      REAL(wp) ::   ztmp, zfwd                 ! Auxiliary variables
1218      LOGICAL  ::   llkbot                     ! Auxiliary variable
1219      !
1220      IF( ln_timing ) CALL timing_start('zdf_osm_vr')
1221      !
1222      zfwd = 1.0_wp
1223      IF( PRESENT(fwd) .AND. ( fwd==.FALSE. ) ) zfwd = -1.0_wp
1224      jktop = 1
1225      IF( PRESENT(ktop) ) jktop = ktop
1226      IF( PRESENT(knlev) ) THEN
1227         jkmax = 0
1228         DO_2D( 0, 0, 0, 0 )
1229            IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj)
1230         END_2D
1231         llkbot = .FALSE.
1232      ELSE
1233         jkmax = jpk
1234         llkbot = .TRUE.
1235      END IF
1236      DO_3D( 0, 0, 0, 0, jktop, jkmax )
1237         IF ( llkbot .OR. knlev(ji,jj) >= jk ) THEN
1238            ztmp         = pu(ji,jj,jk)
1239            pu(ji,jj,jk) = pu(ji,jj,jk) * scos_wind(ji,jj) + zfwd * pv(ji,jj,jk) * ssin_wind(ji,jj)
1240            pv(ji,jj,jk) = pv(ji,jj,jk) * scos_wind(ji,jj) - zfwd * ztmp         * ssin_wind(ji,jj)
1241         END IF
1242      END_3D
1243      !
1244      IF( ln_timing ) CALL timing_stop('zdf_osm_vr')
1245      !
1246   END SUBROUTINE zdf_osm_velocity_rotation_3d
1247
1248   SUBROUTINE zdf_osm_osbl_state( Kmm, pwb_ent, pwb_min, pshear, phbl,   &
1249      &                           phml, pdh )
1250      !!---------------------------------------------------------------------
1251      !!                 ***  ROUTINE zdf_osm_osbl_state  ***
1252      !!
1253      !! ** Purpose : Determines the state of the OSBL, stable/unstable,
1254      !!              shear/ noshear. Also determines shear production,
1255      !!              entrainment buoyancy flux and interfacial Richardson
1256      !!              number
1257      !!
1258      !! ** Method  :
1259      !!
1260      !!----------------------------------------------------------------------
1261      INTEGER,                     INTENT(in   ) ::   Kmm       ! Ocean time-level index
1262      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base
1263      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_min   !    of well-mixed layer
1264      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline
1265      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl      ! BL depth
1266      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phml      ! ML depth
1267      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdh       ! Pycnocline depth
1268      !
1269      ! Local variables
1270      INTEGER :: jj, ji   ! Loop indices
1271      !
1272      REAL(wp), DIMENSION(A2D(0)) ::   zekman
1273      REAL(wp), DIMENSION(A2D(0)) ::   zri_p, zri_b   ! Richardson numbers
1274      REAL(wp)                    ::   zshear_u, zshear_v, zwb_shr
1275      REAL(wp)                    ::   zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
1276      !
1277      REAL(wp), PARAMETER ::   pp_a_shr         = 0.4_wp,  pp_b_shr    = 6.5_wp,  pp_a_wb_s = 0.8_wp
1278      REAL(wp), PARAMETER ::   pp_alpha_c       = 0.2_wp,  pp_alpha_lc = 0.03_wp
1279      REAL(wp), PARAMETER ::   pp_alpha_ls      = 0.06_wp, pp_alpha_s  = 0.15_wp
1280      REAL(wp), PARAMETER ::   pp_ri_p_thresh   = 27.0_wp
1281      REAL(wp), PARAMETER ::   pp_ri_c          = 0.25_wp
1282      REAL(wp), PARAMETER ::   pp_ek            = 4.0_wp
1283      REAL(wp), PARAMETER ::   pp_large         = -1e10_wp
1284      !
1285      IF( ln_timing ) CALL timing_start('zdf_osm_os')
1286      !
1287      ! Initialise arrays
1288      l_conv(:,:)  = .FALSE.
1289      l_shear(:,:) = .FALSE.
1290      n_ddh(:,:)   = 1
1291      ! Initialise INTENT(  out) arrays
1292      pwb_ent(:,:) = pp_large
1293      pwb_min(:,:) = pp_large
1294      !
1295      ! Determins stability and set flag l_conv
1296      DO_2D( 0, 0, 0, 0 )
1297         IF ( shol(ji,jj) < 0._wp ) THEN
1298            l_conv(ji,jj) = .TRUE.
1299         ELSE
1300            l_conv(ji,jj) = .FALSE.
1301         ENDIF
1302      END_2D
1303      !
1304      pshear(:,:) = 0.0_wp
1305      zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D(0)) ) * phbl(:,:) / MAX( sustar(A2D(0)), 1.e-8 ) )
1306      !
1307      DO_2D( 0, 0, 0, 0 )
1308         IF ( l_conv(ji,jj) ) THEN
1309            IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN
1310               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,       &
1311                  &                                                          1e-8_wp ) ) * ( phbl(ji,jj) / pdh(ji,jj) ) *   &
1312                  &                  ( svstr(ji,jj) / MAX( sustar(ji,jj), 1e-6_wp ) )**2 /                                  &
1313                  &                  MAX( zekman(ji,jj), 1.0e-6_wp ), 5.0_wp )
1314               IF ( ff_t(ji,jj) >= 0.0_wp ) THEN   ! Northern hemisphere
1315                  zri_b(ji,jj) = av_db_ml(ji,jj) * pdh(ji,jj) / ( MAX( av_du_ml(ji,jj), 1e-5_wp )**2 +   &
1316                     &                                          MAX( -1.0_wp * av_dv_ml(ji,jj), 1e-5_wp)**2 )
1317               ELSE                                ! Southern hemisphere
1318                  zri_b(ji,jj) = av_db_ml(ji,jj) * pdh(ji,jj) / ( MAX( av_du_ml(ji,jj), 1e-5_wp )**2 +   &
1319                     &                                          MAX(           av_dv_ml(ji,jj), 1e-5_wp)**2 )
1320               END IF
1321               pshear(ji,jj) = pp_a_shr * zekman(ji,jj) *                                                   &
1322                  &            ( MAX( sustar(ji,jj)**2 * av_du_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) +          &
1323                  &              pp_b_shr * MAX( -1.0_wp * ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) *   &
1324                  &                            av_dv_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) )
1325               ! Stability dependence
1326               pshear(ji,jj) = pshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - pp_ri_c ) / pp_ri_c ) )
1327               !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1328               ! Test ensures n_ddh=0 is not selected. Change to zri_p<27 when  !
1329               ! full code available                                          !
1330               !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1331               IF ( pshear(ji,jj) > 1e-10 ) THEN
1332                  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
1333                     ! Growing shear layer
1334                     n_ddh(ji,jj) = 0
1335                     l_shear(ji,jj) = .TRUE.
1336                  ELSE
1337                     n_ddh(ji,jj) = 1
1338                     !             IF ( zri_b <= 1.5 .and. pshear(ji,jj) > 0._wp ) THEN
1339                     ! Shear production large enough to determine layer charcteristics, but can't maintain a shear layer
1340                     l_shear(ji,jj) = .TRUE.
1341                     !             ELSE
1342                  END IF
1343               ELSE
1344                  n_ddh(ji,jj) = 2
1345                  l_shear(ji,jj) = .FALSE.
1346               END IF
1347               ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline
1348               !               pshear(ji,jj) = 0.5 * pshear(ji,jj)
1349               !               l_shear(ji,jj) = .FALSE.
1350               !            ENDIF
1351            ELSE   ! av_db_bl test, note pshear set to zero
1352               n_ddh(ji,jj) = 2
1353               l_shear(ji,jj) = .FALSE.
1354            ENDIF
1355         ENDIF
1356      END_2D
1357      !
1358      ! Calculate entrainment buoyancy flux due to surface fluxes.
1359      DO_2D( 0, 0, 0, 0 )
1360         IF ( l_conv(ji,jj) ) THEN
1361            zwcor        = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln
1362            zrf_conv     = TANH( ( swstrc(ji,jj) / zwcor )**0.69_wp )
1363            zrf_shear    = TANH( ( sustar(ji,jj) / zwcor )**0.69_wp )
1364            zrf_langmuir = TANH( ( swstrl(ji,jj) / zwcor )**0.69_wp )
1365            IF ( nn_osm_SD_reduce > 0 ) THEN
1366               ! Effective Stokes drift already reduced from surface value
1367               zr_stokes = 1.0_wp
1368            ELSE
1369               ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
1370               ! requires further reduction where BL is deep
1371               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) ) )
1372            END IF
1373            pwb_ent(ji,jj) = -2.0_wp * pp_alpha_c * zrf_conv * swbav(ji,jj) -                                          &
1374               &             pp_alpha_s * zrf_shear * sustar(ji,jj)**3 / phml(ji,jj) +                                 &
1375               &             zr_stokes * ( pp_alpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 -   &
1376               &                           zrf_langmuir * pp_alpha_lc * swstrl(ji,jj)**3 ) / phml(ji,jj)
1377         ENDIF
1378      END_2D
1379      !
1380      DO_2D( 0, 0, 0, 0 )
1381         IF ( l_shear(ji,jj) ) THEN
1382            IF ( l_conv(ji,jj) ) THEN
1383               ! Unstable OSBL
1384               zwb_shr = -1.0_wp * pp_a_wb_s * zri_b(ji,jj) * pshear(ji,jj)
1385               IF ( n_ddh(ji,jj) == 0 ) THEN
1386                  ! Developing shear layer, additional shear production possible.
1387
1388                  !              pshear_u = MAX( zustar(ji,jj)**2 * MAX( av_du_ml(ji,jj), 0._wp ) /  phbl(ji,jj), 0._wp )
1389                  !              pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1.d0 )**2 )
1390                  !              pshear(ji,jj) = MIN( pshear(ji,jj), pshear_u )
1391
1392                  !              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 )
1393                  !              zwb_shr = MAX( zwb_shr, -0.25 * pshear_u )
1394               ENDIF
1395               pwb_ent(ji,jj) = pwb_ent(ji,jj) + zwb_shr
1396               !           pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * zwb0(ji,jj)
1397            ELSE   ! IF ( l_conv ) THEN - ENDIF
1398               ! Stable OSBL  - shear production not coded for first attempt.
1399            ENDIF   ! l_conv
1400         END IF   ! l_shear
1401         IF ( l_conv(ji,jj) ) THEN
1402            ! Unstable OSBL
1403            pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * 2.0_wp * swbav(ji,jj)
1404         END IF  ! l_conv
1405      END_2D
1406      !
1407      IF( ln_timing ) CALL timing_stop('zdf_osm_os')
1408      !
1409   END SUBROUTINE zdf_osm_osbl_state
1410
1411   SUBROUTINE zdf_osm_external_gradients( Kmm, kbase, pdtdz, pdsdz, pdbdz )
1412      !!---------------------------------------------------------------------
1413      !!                   ***  ROUTINE zdf_osm_external_gradients  ***
1414      !!
1415      !! ** Purpose : Calculates the gradients below the OSBL
1416      !!
1417      !! ** Method  : Uses nbld and ibld_ext to determine levels to calculate the gradient.
1418      !!
1419      !!----------------------------------------------------------------------   
1420      INTEGER,                     INTENT(in   ) ::   Kmm                   ! Ocean time-level index
1421      INTEGER,  DIMENSION(A2D(0)), INTENT(in   ) ::   kbase                 ! OSBL base layer index
1422      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdtdz, pdsdz, pdbdz   ! External gradients of temperature, salinity and buoyancy
1423      !
1424      ! Local variables
1425      INTEGER  ::   ji, jj, jkb, jkb1
1426      REAL(wp) ::   zthermal, zbeta
1427      !
1428      REAL(wp), PARAMETER ::   pp_large = -1e10_wp
1429      !
1430      IF( ln_timing ) CALL timing_start('zdf_osm_eg')
1431      !
1432      pdtdz(:,:) = pp_large
1433      pdsdz(:,:) = pp_large
1434      pdbdz(:,:) = pp_large
1435      !
1436      DO_2D( 0, 0, 0, 0 )
1437         IF ( kbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
1438            zthermal = rab_n(ji,jj,1,jp_tem)   ! Ideally use nbld not 1??
1439            zbeta    = rab_n(ji,jj,1,jp_sal)
1440            jkb = kbase(ji,jj)
1441            jkb1 = MIN( jkb + 1, mbkt(ji,jj) )
1442            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)
1443            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)
1444            pdbdz(ji,jj) = grav * zthermal * pdtdz(ji,jj) - grav * zbeta * pdsdz(ji,jj)
1445         ELSE
1446            pdtdz(ji,jj) = 0.0_wp
1447            pdsdz(ji,jj) = 0.0_wp
1448            pdbdz(ji,jj) = 0.0_wp
1449         END IF
1450      END_2D
1451      !
1452      IF( ln_timing ) CALL timing_stop('zdf_osm_eg')
1453      !
1454   END SUBROUTINE zdf_osm_external_gradients
1455
1456   SUBROUTINE zdf_osm_calculate_dhdt( pdhdt, phbl, pdh, pwb_ent, pwb_min,   &
1457      &                               pdbdz_bl_ext, pwb_fk_b, pwb_fk, pvel_mle )
1458      !!---------------------------------------------------------------------
1459      !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
1460      !!
1461      !! ** Purpose : Calculates the rate at which hbl changes.
1462      !!
1463      !! ** Method  :
1464      !!
1465      !!----------------------------------------------------------------------
1466      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdhdt          ! Rate of change of hbl
1467      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl           ! BL depth
1468      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdh            ! Pycnocline depth
1469      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
1470      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_min
1471      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
1472      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL
1473      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux
1474      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK
1475      !
1476      ! Local variables
1477      INTEGER  ::   jj, ji
1478      REAL(wp) ::   zgamma_b_nd, zgamma_dh_nd, zpert, zpsi, zari
1479      REAL(wp) ::   zvel_max, zddhdt
1480      !
1481      REAL(wp), PARAMETER ::   pp_alpha_b = 0.3_wp
1482      REAL(wp), PARAMETER ::   pp_ddh     = 2.5_wp, pp_ddh_2 = 3.5_wp   ! Also in pycnocline_depth
1483      REAL(wp), PARAMETER ::   pp_large   = -1e10_wp
1484      !
1485      IF( ln_timing ) CALL timing_start('zdf_osm_cd')
1486      !
1487      pdhdt(:,:)    = pp_large
1488      pwb_fk_b(:,:) = pp_large
1489      !
1490      DO_2D( 0, 0, 0, 0 )
1491         !
1492         IF ( l_shear(ji,jj) ) THEN
1493            !
1494            IF ( l_conv(ji,jj) ) THEN   ! Convective
1495               !
1496               IF ( ln_osm_mle ) THEN
1497                  IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN   ! Fox-Kemper buoyancy flux average over OSBL
1498                     pwb_fk_b(ji,jj) = pwb_fk(ji,jj) * ( 1.0_wp + hmle(ji,jj) / ( 6.0_wp * hbl(ji,jj) ) *   &
1499                        &                                         ( -1.0_wp + ( 1.0_wp - 2.0_wp * hbl(ji,jj) / hmle(ji,jj) )**3 ) )
1500                  ELSE
1501                     pwb_fk_b(ji,jj) = 0.5_wp * pwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
1502                  ENDIF
1503                  zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1504                  IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN   ! OSBL is deepening,
1505                     !                                                                 !    entrainment > restratification
1506                     IF ( av_db_bl(ji,jj) > 1e-15_wp ) THEN
1507                        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 ) )
1508                        zpsi = ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   &
1509                           &   ( 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)
1510                        zpsi = zpsi + 1.75_wp * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   &
1511                           &   ( 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 )
1512                        zpsi = pp_alpha_b * MAX( zpsi, 0.0_wp )
1513                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /   &
1514                           &                      ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) +   &
1515                           &            zpsi / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
1516                        IF ( n_ddh(ji,jj) == 1 ) THEN
1517                           IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5_wp ) THEN
1518                              zari = MIN( 1.5_wp * av_db_bl(ji,jj) /                                                   &
1519                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                     &
1520                                 &                               av_db_bl(ji,jj)**2 / MAX( 4.5_wp * svstr(ji,jj)**2,   &
1521                                 &                                                       1e-12_wp ) ) ), 0.2_wp )
1522                           ELSE
1523                              zari = MIN( 1.5_wp * av_db_bl(ji,jj) /                                                    &
1524                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                      &
1525                                 &                               av_db_bl(ji,jj)**2 / MAX( 4.5_wp * swstrc(ji,jj)**2,   &
1526                                 &                                                       1e-12_wp ) ) ), 0.2_wp )
1527                           ENDIF
1528                           ! Relaxation to dh_ref = zari * hbl
1529                           zddhdt = -1.0_wp * pp_ddh_2 * ( 1.0_wp - pdh(ji,jj) / ( zari * phbl(ji,jj) ) ) * pwb_ent(ji,jj) /   &
1530                              &     ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
1531                        ELSE IF ( n_ddh(ji,jj) == 0 ) THEN   ! Growing shear layer
1532                           zddhdt = -1.0_wp * pp_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) /   &
1533                              &     ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
1534                           zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX( sustar(ji,jj), 1e-8_wp ) ) * zddhdt
1535                        ELSE
1536                           zddhdt = 0.0_wp
1537                        ENDIF   ! n_ddh
1538                        pdhdt(ji,jj) = pdhdt(ji,jj) + pp_alpha_b * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   &
1539                           &                            av_db_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) /   &
1540                           &                            ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
1541                     ELSE   ! av_db_bl >0
1542                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1e-15_wp )
1543                     ENDIF
1544                  ELSE   ! pwb_min + 2*pwb_fk_b < 0
1545                     ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
1546                     pdhdt(ji,jj) = -1.0_wp * MIN( pvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp )
1547                  ENDIF
1548               ELSE   ! Fox-Kemper not used.
1549                  zvel_max = -1.0_wp * ( 1.0_wp + 1.0_wp * ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird *     &
1550                     &                                                         rn_Dt / hbl(ji,jj) ) * pwb_ent(ji,jj) /   &
1551                     &       MAX( ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln )
1552                  pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
1553                  ! added ajgn 23 July as temporay fix
1554               ENDIF   ! ln_osm_mle
1555               !
1556            ELSE   ! l_conv - Stable
1557               !
1558               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)
1559               IF ( pdhdt(ji,jj) < 0.0_wp ) THEN   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
1560                  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)
1561               ELSE
1562                  zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), av_db_bl(ji,jj) )
1563               ENDIF
1564               pdhdt(ji,jj) = 2.0_wp * pdhdt(ji,jj) / MAX( zpert, epsln )
1565               pdhdt(ji,jj) = MAX( pdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp )
1566               !
1567            ENDIF   ! l_conv
1568            !
1569         ELSE   ! l_shear
1570            !
1571            IF ( l_conv(ji,jj) ) THEN   ! Convective
1572               !
1573               IF ( ln_osm_mle ) THEN
1574                  IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN   ! Fox-Kemper buoyancy flux average over OSBL
1575                     pwb_fk_b(ji,jj) = pwb_fk(ji,jj) *                       &
1576                        ( 1.0_wp + hmle(ji,jj) / ( 6.0_wp * hbl(ji,jj) ) *   &
1577                        &          ( -1.0_wp + ( 1.0_wp - 2.0_wp * hbl(ji,jj) / hmle(ji,jj))**3) )
1578                  ELSE
1579                     pwb_fk_b(ji,jj) = 0.5_wp * pwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
1580                  ENDIF
1581                  zvel_max = ( swstrl(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1582                  IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN   ! OSBL is deepening,
1583                     !                                                                 !    entrainment > restratification
1584                     IF ( av_db_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN
1585                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /   &
1586                           &            ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
1587                     ELSE
1588                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) / MAX( zvel_max, 1e-15_wp )
1589                     ENDIF
1590                  ELSE   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
1591                     pdhdt(ji,jj) = -1.0_wp * MIN( pvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp )
1592                  ENDIF
1593               ELSE   ! Fox-Kemper not used
1594                  zvel_max = -1.0_wp * pwb_ent(ji,jj) / MAX( ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln )
1595                  pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
1596                  ! added ajgn 23 July as temporay fix
1597               ENDIF  ! ln_osm_mle
1598               !
1599            ELSE                        ! Stable
1600               !
1601               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)
1602               IF ( pdhdt(ji,jj) < 0.0_wp ) THEN
1603                  ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
1604                  zpert = 2.0_wp * svstr(ji,jj)**2 / hbl(ji,jj)
1605               ELSE
1606                  zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), av_db_bl(ji,jj) )
1607               ENDIF
1608               pdhdt(ji,jj) = 2.0_wp * pdhdt(ji,jj) / MAX(zpert, epsln)
1609               pdhdt(ji,jj) = MAX( pdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp )
1610               !
1611            ENDIF  ! l_conv
1612            !
1613         ENDIF ! l_shear
1614         !
1615      END_2D
1616      !
1617      IF( ln_timing ) CALL timing_stop('zdf_osm_cd')
1618      !
1619   END SUBROUTINE zdf_osm_calculate_dhdt
1620
1621   SUBROUTINE zdf_osm_timestep_hbl( Kmm, pdhdt, phbl, phbl_t, pwb_ent,   &
1622      &                             pwb_fk_b )
1623      !!---------------------------------------------------------------------
1624      !!                ***  ROUTINE zdf_osm_timestep_hbl  ***
1625      !!
1626      !! ** Purpose : Increments hbl.
1627      !!
1628      !! ** Method  : If the change in hbl exceeds one model level the change is
1629      !!              is calculated by moving down the grid, changing the
1630      !!              buoyancy jump. This is to ensure that the change in hbl
1631      !!              does not overshoot a stable layer.
1632      !!
1633      !!----------------------------------------------------------------------
1634      INTEGER,                     INTENT(in   ) ::   Kmm        ! Ocean time-level index
1635      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdhdt      ! Rates of change of hbl
1636      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phbl       ! BL depth
1637      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl_t     ! BL depth
1638      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux
1639      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL
1640      !
1641      ! Local variables
1642      INTEGER  ::   jk, jj, ji, jm
1643      REAL(wp) ::   zhbl_s, zvel_max, zdb
1644      REAL(wp) ::   zthermal, zbeta
1645      !
1646      IF( ln_timing ) CALL timing_start('zdf_osm_th')
1647      !
1648      DO_2D( 0, 0, 0, 0 )
1649         IF ( nbld(ji,jj) - nmld(ji,jj) > 1 ) THEN
1650            !
1651            ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
1652            !
1653            zhbl_s   = hbl(ji,jj)
1654            jm       = nmld(ji,jj)
1655            zthermal = rab_n(ji,jj,1,jp_tem)
1656            zbeta    = rab_n(ji,jj,1,jp_sal)
1657            !
1658            IF ( l_conv(ji,jj) ) THEN   ! Unstable
1659               !
1660               IF( ln_osm_mle ) THEN
1661                  zvel_max = ( swstrl(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1662               ELSE
1663                  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 /   &
1664                     &                                     hbl(ji,jj) ) * pwb_ent(ji,jj) /                                     &
1665                     &       ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird
1666               ENDIF
1667               DO jk = nmld(ji,jj), nbld(ji,jj)
1668                  zdb = MAX( grav * ( zthermal * ( av_t_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) -   &
1669                     &                zbeta    * ( av_s_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) + zvel_max
1670                  !
1671                  IF ( ln_osm_mle ) THEN
1672                     zhbl_s = zhbl_s + MIN( rn_Dt * ( ( -1.0_wp * pwb_ent(ji,jj) - 2.0_wp * pwb_fk_b(ji,jj) ) / zdb ) /   &
1673                        &                   REAL( nbld(ji,jj) - nmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) )
1674                  ELSE
1675                     zhbl_s = zhbl_s + MIN( rn_Dt * ( -1.0_wp * pwb_ent(ji,jj) / zdb ) /   &
1676                        &                   REAL( nbld(ji,jj) - nmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) )
1677                  ENDIF
1678                  !                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
1679                  IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN
1680                     zhbl_s = MIN( zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm ) - depth_tol )
1681                     l_pyc(ji,jj) = .FALSE.
1682                  ENDIF
1683                  IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1
1684               END DO
1685               hbl(ji,jj)  = zhbl_s
1686               nbld(ji,jj) = jm
1687            ELSE   ! Stable
1688               DO jk = nmld(ji,jj), nbld(ji,jj)
1689                  zdb = MAX(  grav * ( zthermal * ( av_t_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) -               &
1690                     &                 zbeta    * ( av_s_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) +   &
1691                     &  2.0 * svstr(ji,jj)**2 / zhbl_s
1692                  !
1693                  ! Alan is thuis right? I have simply changed hbli to hbl
1694                  shol(ji,jj)  = -1.0_wp * zhbl_s / ( ( svstr(ji,jj)**3 + epsln ) / swbav(ji,jj) )
1695                  pdhdt(ji,jj) = -1.0_wp * ( swbav(ji,jj) - 0.04_wp / 2.0_wp * swstrl(ji,jj)**3 / zhbl_s -   &
1696                     &                       0.15_wp / 2.0_wp * ( 1.0_wp - EXP( -1.5_wp * sla(ji,jj) ) ) *   &
1697                     &                                 sustar(ji,jj)**3 / zhbl_s ) *                         &
1698                     &           ( 0.725_wp + 0.225_wp * EXP( -7.5_wp * shol(ji,jj) ) )
1699                  pdhdt(ji,jj) = pdhdt(ji,jj) + swbav(ji,jj)
1700                  zhbl_s = zhbl_s + MIN( pdhdt(ji,jj) / zdb * rn_Dt / REAL( nbld(ji,jj) - nmld(ji,jj), KIND=wp ),   &
1701                     &                   e3w(ji,jj,jm,Kmm) )
1702                 
1703                  !                    zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
1704                  IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN
1705                     zhbl_s      = MIN( zhbl_s,  gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - depth_tol )
1706                     l_pyc(ji,jj) = .FALSE.
1707                  ENDIF
1708                  IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1
1709               END DO
1710            ENDIF   ! IF ( l_conv )
1711            hbl(ji,jj)  = MAX( zhbl_s, gdepw(ji,jj,4,Kmm) )
1712            nbld(ji,jj) = MAX( jm, 4 )
1713         ELSE
1714            ! change zero or one model level.
1715            hbl(ji,jj) = MAX( phbl_t(ji,jj), gdepw(ji,jj,4,Kmm) )
1716         ENDIF
1717         phbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm)
1718      END_2D
1719      !
1720      IF( ln_timing ) CALL timing_stop('zdf_osm_th')
1721      !
1722   END SUBROUTINE zdf_osm_timestep_hbl
1723
1724   SUBROUTINE zdf_osm_pycnocline_thickness( Kmm, pdh, phml, pdhdt, phbl,   &
1725      &                                     pwb_ent, pdbdz_bl_ext, pwb_fk_b )
1726      !!---------------------------------------------------------------------
1727      !!            ***  ROUTINE zdf_osm_pycnocline_thickness  ***
1728      !!
1729      !! ** Purpose : Calculates thickness of the pycnocline
1730      !!
1731      !! ** Method  : The thickness is calculated from a prognostic equation
1732      !!              that relaxes the pycnocine thickness to a diagnostic
1733      !!              value. The time change is calculated assuming the
1734      !!              thickness relaxes exponentially. This is done to deal
1735      !!              with large timesteps.
1736      !!
1737      !!----------------------------------------------------------------------
1738      INTEGER,                     INTENT(in   ) ::   Kmm            ! Ocean time-level index
1739      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdh            ! Pycnocline thickness
1740      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phml           ! ML depth
1741      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency
1742      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl           ! BL depth
1743      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
1744      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
1745      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL
1746
1747      !
1748      ! Local variables
1749      INTEGER  ::   jj, ji
1750      INTEGER  ::   inhml
1751      REAL(wp) ::   zari, ztau, zdh_ref, zddhdt, zvel_max
1752      REAL(wp) ::   ztmp   ! Auxiliary variable
1753      !
1754      REAL, PARAMETER ::   pp_ddh = 2.5_wp, pp_ddh_2 = 3.5_wp   ! Also in pycnocline_depth
1755      !
1756      IF( ln_timing ) CALL timing_start('zdf_osm_pt')
1757      !
1758      DO_2D( 0, 0, 0, 0 )
1759         !
1760         IF ( l_shear(ji,jj) ) THEN
1761            !
1762            IF ( l_conv(ji,jj) ) THEN
1763               !
1764               IF ( av_db_bl(ji,jj) > 1e-15_wp ) THEN
1765                  IF ( n_ddh(ji,jj) == 0 ) THEN
1766                     zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1767                     ! ddhdt for pycnocline determined in osm_calculate_dhdt
1768                     zddhdt = -1.0_wp * pp_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) /   &
1769                        &     ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15 ) )
1770                     zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX( sustar(ji,jj), 1e-8 ) ) * zddhdt
1771                     ! Maximum limit for how thick the shear layer can grow relative to the thickness of the boundary layer
1772                     dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_Dt, 0.625_wp * hbl(ji,jj) )
1773                  ELSE   ! Need to recalculate because hbl has been updated
1774                     IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5_wp ) THEN
1775                        ztmp = svstr(ji,jj)
1776                     ELSE
1777                        ztmp = swstrc(ji,jj)
1778                     END IF
1779                     zari = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +        &
1780                        &                                                   av_db_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2,   &
1781                        &                                                                           1e-12_wp ) ) ), 0.2_wp )
1782                     ztau = MAX( av_db_bl(ji,jj) * ( zari * hbl(ji,jj) ) /   &
1783                        &        ( pp_ddh_2 * MAX( -1.0_wp * pwb_ent(ji,jj), 1e-12_wp ) ), 2.0_wp * rn_Dt )
1784                     dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) +   &
1785                        &        zari * phbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) )
1786                     IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * phbl(ji,jj)
1787                  END IF
1788               ELSE
1789                  ztau = MAX( MAX( hbl(ji,jj) / ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln), 2.0_wp * rn_Dt )
1790                  dh(ji,jj) = dh(ji,jj) * EXP( -1.0_wp * rn_Dt / ztau ) +   &
1791                     &        0.2_wp * phbl(ji,jj) * ( 1.0_wp - EXP( -1.0_wp * rn_Dt / ztau ) )
1792                  IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2_wp * hbl(ji,jj)
1793               END IF
1794               !
1795            ELSE   ! l_conv
1796               ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL
1797               ztau = hbl(ji,jj) / MAX(svstr(ji,jj), epsln)
1798               IF ( pdhdt(ji,jj) >= 0.0_wp ) THEN   ! Probably shouldn't include wm here
1799                  ! Boundary layer deepening
1800                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN
1801                     ! Pycnocline thickness set by stratification - use same relationship as for neutral conditions
1802                     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 )
1803                     zdh_ref = MIN( zari, 0.2_wp ) * hbl(ji,jj)
1804                  ELSE
1805                     zdh_ref = 0.2_wp * hbl(ji,jj)
1806                  ENDIF
1807               ELSE   ! IF(dhdt < 0)
1808                  zdh_ref = 0.2_wp * hbl(ji,jj)
1809               ENDIF   ! IF (dhdt >= 0)
1810               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 ) )
1811               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
1812               !                                                                                !    rapid collapse
1813            ENDIF
1814            !
1815         ELSE   ! l_shear = .FALSE., calculate ddhdt here
1816            !
1817            IF ( l_conv(ji,jj) ) THEN
1818               !
1819               IF( ln_osm_mle ) THEN
1820                  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
1821                     !                                                                 !    ln_osm_mle=F
1822                     IF ( av_db_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN
1823                        IF ( ( swstrc(ji,jj) / MAX( svstr(ji,jj), epsln) )**3 <= 0.5_wp ) THEN   ! Near neutral stability
1824                           ztmp = svstr(ji,jj)
1825                        ELSE   ! Unstable
1826                           ztmp = swstrc(ji,jj)
1827                        END IF
1828                        zari = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) *                           &
1829                           &                                   ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp) +   &
1830                           &                                     av_db_bl(ji,jj)**2 / MAX(4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp )
1831                     ELSE
1832                        zari = 0.2_wp
1833                     END IF
1834                  ELSE
1835                     zari = 0.2_wp
1836                  END IF
1837                  ztau    = 0.2_wp * hbl(ji,jj) / MAX( epsln, ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird )
1838                  zdh_ref = zari * hbl(ji,jj)
1839               ELSE   ! ln_osm_mle
1840                  IF ( av_db_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN
1841                     IF ( ( swstrc(ji,jj) / MAX( svstr(ji,jj), epsln ) )**3 <= 0.5_wp ) THEN   ! Near neutral stability
1842                        ztmp = svstr(ji,jj)
1843                     ELSE   ! Unstable
1844                        ztmp = swstrc(ji,jj)
1845                     END IF
1846                     zari    = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) *                            &
1847                        &                                      ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +   &
1848                        &                                        av_db_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp )
1849                  ELSE
1850                     zari    = 0.2_wp
1851                  END IF
1852                  ztau    = hbl(ji,jj) / MAX( epsln, ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird )
1853                  zdh_ref = zari * hbl(ji,jj)
1854               END IF   ! ln_osm_mle
1855               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 ) )
1856               !               IF ( pdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
1857               IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
1858               ! Alan: this hml is never defined or used
1859            ELSE   ! IF (l_conv)
1860               !
1861               ztau = hbl(ji,jj) / MAX( svstr(ji,jj), epsln )
1862               IF ( pdhdt(ji,jj) >= 0.0_wp ) THEN   ! Probably shouldn't include wm here
1863                  ! Boundary layer deepening
1864                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN
1865                     ! Pycnocline thickness set by stratification - use same relationship as for neutral conditions.
1866                     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 )
1867                     zdh_ref = MIN( zari, 0.2_wp ) * hbl(ji,jj)
1868                  ELSE
1869                     zdh_ref = 0.2_wp * hbl(ji,jj)
1870                  END IF
1871               ELSE   ! IF(dhdt < 0)
1872                  zdh_ref = 0.2_wp * hbl(ji,jj)
1873               END IF   ! IF (dhdt >= 0)
1874               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 ) )
1875               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
1876               !                                                                                !    rapid collapse
1877            END IF   ! IF (l_conv)
1878            !
1879         END IF   ! l_shear
1880         !
1881         hml(ji,jj)  = hbl(ji,jj) - dh(ji,jj)
1882         inhml       = MAX( INT( dh(ji,jj) / MAX( e3t(ji,jj,nbld(ji,jj)-1,Kmm), 1e-3_wp ) ), 1 )
1883         nmld(ji,jj) = MAX( nbld(ji,jj) - inhml, 3 )
1884         phml(ji,jj) = gdepw(ji,jj,nmld(ji,jj),Kmm)
1885         pdh(ji,jj)  = phbl(ji,jj) - phml(ji,jj)
1886         !
1887      END_2D
1888      !
1889      IF( ln_timing ) CALL timing_stop('zdf_osm_pt')
1890      !
1891   END SUBROUTINE zdf_osm_pycnocline_thickness
1892
1893   SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( Kmm, kp_ext, pdbdz, palpha, pdh,   &
1894      &                                             phbl, pdbdz_bl_ext, phml, pdhdt )
1895      !!---------------------------------------------------------------------
1896      !!       ***  ROUTINE zdf_osm_pycnocline_buoyancy_profiles  ***
1897      !!
1898      !! ** Purpose : calculate pycnocline buoyancy profiles
1899      !!
1900      !! ** Method  :
1901      !!
1902      !!----------------------------------------------------------------------
1903      INTEGER,                          INTENT(in   ) ::   Kmm            ! Ocean time-level index
1904      INTEGER,  DIMENSION(A2D(0)),      INTENT(in   ) ::   kp_ext         ! External-level offsets
1905      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(  out) ::   pdbdz          ! Gradients in the pycnocline
1906      REAL(wp), DIMENSION(A2D(0)),      INTENT(  out) ::   palpha
1907      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline thickness
1908      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth
1909      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
1910      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth
1911      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! Rates of change of hbl
1912      !
1913      ! Local variables
1914      INTEGER  ::   jk, jj, ji
1915      REAL(wp) ::   zbgrad
1916      REAL(wp) ::   zgamma_b_nd, znd
1917      REAL(wp) ::   zzeta_m
1918      REAL(wp) ::   ztmp   ! Auxiliary variable
1919      !
1920      REAL(wp), PARAMETER ::   pp_gamma_b = 2.25_wp
1921      REAL(wp), PARAMETER ::   pp_large   = -1e10_wp
1922      !
1923      IF( ln_timing ) CALL timing_start('zdf_osm_pscp')
1924      !
1925      pdbdz(:,:,:) = pp_large
1926      palpha(:,:)  = pp_large
1927      !
1928      DO_2D( 0, 0, 0, 0 )
1929         !
1930         IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN
1931            !
1932            IF ( l_conv(ji,jj) ) THEN   ! Convective conditions
1933               !
1934               IF ( l_pyc(ji,jj) ) THEN
1935                  !
1936                  zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )
1937                  palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / pp_gamma_b ) ) *   &
1938                     &                                pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / av_db_ml(ji,jj) ) /                &
1939                     &            ( 0.723_wp + SQRT( 3.14159_wp / pp_gamma_b ) )
1940                  palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp )
1941                  ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln )
1942                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1943                  ! Commented lines in this section are not needed in new code, once tested !
1944                  ! can be removed                                                          !
1945                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1946                  ! ztgrad = zalpha * av_dt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
1947                  ! zsgrad = zalpha * av_ds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
1948                  zbgrad = palpha(ji,jj) * av_db_ml(ji,jj) * ztmp + pdbdz_bl_ext(ji,jj)
1949                  zgamma_b_nd = pdbdz_bl_ext(ji,jj) * pdh(ji,jj) / MAX( av_db_ml(ji,jj), epsln )
1950                  DO jk = 2, nbld(ji,jj)
1951                     znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) * ztmp
1952                     IF ( znd <= zzeta_m ) THEN
1953                        ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * av_dt_ml(ji,jj) * ztmp * &
1954                        ! &        EXP( -6.0 * ( znd -zzeta_m )**2 )
1955                        ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * av_ds_ml(ji,jj) * ztmp * &
1956                        ! & EXP( -6.0 * ( znd -zzeta_m )**2 )
1957                        pdbdz(ji,jj,jk) = pdbdz_bl_ext(ji,jj) + palpha(ji,jj) * av_db_ml(ji,jj) * ztmp * &
1958                           & EXP( -6.0_wp * ( znd -zzeta_m )**2 )
1959                     ELSE
1960                        ! zdtdz(ji,jj,jk) =  ztgrad * EXP( -pp_gamma_b * ( znd - zzeta_m )**2 )
1961                        ! zdsdz(ji,jj,jk) =  zsgrad * EXP( -pp_gamma_b * ( znd - zzeta_m )**2 )
1962                        pdbdz(ji,jj,jk) =  zbgrad * EXP( -1.0_wp * pp_gamma_b * ( znd - zzeta_m )**2 )
1963                     END IF
1964                  END DO
1965               END IF   ! If no pycnocline pycnocline gradients set to zero
1966               !
1967            ELSE   ! Stable conditions
1968               ! If pycnocline profile only defined when depth steady of increasing.
1969               IF ( pdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady.
1970                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN
1971                     IF ( shol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline
1972                        ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln )
1973                        zbgrad = av_db_bl(ji,jj) * ztmp
1974                        DO jk = 2, nbld(ji,jj)
1975                           znd = gdepw(ji,jj,jk,Kmm) * ztmp
1976                           pdbdz(ji,jj,jk) = zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 )
1977                        END DO
1978                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
1979                        ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln )
1980                        zbgrad = av_db_bl(ji,jj) * ztmp
1981                        DO jk = 2, nbld(ji,jj)
1982                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp
1983                           pdbdz(ji,jj,jk) = zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 )
1984                        END DO
1985                     END IF   ! IF (shol >=0.5)
1986                  END IF      ! IF (av_db_bl> 0.)
1987               END IF         ! IF (pdhdt >= 0) pdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
1988               !
1989            END IF            ! IF (l_conv)
1990            !
1991         END IF   ! IF ( nbld(ji,jj) < mbkt(ji,jj) )
1992         !
1993      END_2D
1994      !
1995      IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles
1996         IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) )
1997      END IF
1998      !
1999      IF( ln_timing ) CALL timing_stop('zdf_osm_pscp')
2000      !
2001   END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles
2002
2003   SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, pdiffut, pviscos, phbl,   &
2004      &                                      phml, pdh, pdhdt, pshear,           &
2005      &                                      pwb_ent, pwb_min )
2006      !!---------------------------------------------------------------------
2007      !!           ***  ROUTINE zdf_osm_diffusivity_viscosity  ***
2008      !!
2009      !! ** Purpose : Determines the eddy diffusivity and eddy viscosity
2010      !!              profiles in the mixed layer and the pycnocline.
2011      !!
2012      !! ** Method  :
2013      !!
2014      !!----------------------------------------------------------------------
2015      INTEGER,                          INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices
2016      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(inout) ::   pdiffut        ! t-diffusivity
2017      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(inout) ::   pviscos        ! Viscosity
2018      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth
2019      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth
2020      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline depth
2021      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency
2022      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pshear         ! Shear production
2023      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
2024      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pwb_min
2025      !
2026      ! Local variables
2027      INTEGER ::   ji, jj, jk   ! Loop indices
2028      !
2029      ! Scales used to calculate eddy diffusivity and viscosity profiles
2030      REAL(wp), DIMENSION(A2D(0)) ::   zdifml_sc,    zvisml_sc
2031      REAL(wp), DIMENSION(A2D(0)) ::   zdifpyc_n_sc, zdifpyc_s_sc
2032      REAL(wp), DIMENSION(A2D(0)) ::   zvispyc_n_sc, zvispyc_s_sc
2033      REAL(wp), DIMENSION(A2D(0)) ::   zbeta_d_sc,   zbeta_v_sc
2034      REAL(wp), DIMENSION(A2D(0)) ::   zb_coup,      zc_coup_vis,  zc_coup_dif
2035      !
2036      REAL(wp) ::   zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b
2037      REAL(wp) ::   za_cubic, zb_d_cubic, zc_d_cubic, zd_d_cubic,   &   ! Coefficients in cubic polynomial specifying diffusivity
2038         &                    zb_v_cubic, zc_v_cubic, zd_v_cubic        ! and viscosity in pycnocline
2039      REAL(wp) ::   zznd_ml, zznd_pyc, ztmp
2040      REAL(wp) ::   zmsku, zmskv
2041      !
2042      REAL(wp), PARAMETER ::   pp_dif_ml     = 0.8_wp,  pp_vis_ml  = 0.375_wp
2043      REAL(wp), PARAMETER ::   pp_dif_pyc    = 0.15_wp, pp_vis_pyc = 0.142_wp
2044      REAL(wp), PARAMETER ::   pp_vispyc_shr = 0.15_wp
2045      !
2046      IF( ln_timing ) CALL timing_start('zdf_osm_dv')
2047      !
2048      zb_coup(:,:) = 0.0_wp
2049      !
2050      DO_2D( 0, 0, 0, 0 )
2051         IF ( l_conv(ji,jj) ) THEN
2052            !
2053            zvel_sc_pyc = ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 + 4.25_wp * pshear(ji,jj) * phbl(ji,jj) )**pthird
2054            zvel_sc_ml  = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird
2055            zstab_fac   = ( phml(ji,jj) / zvel_sc_ml *   &
2056               &            ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP(-3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )**1.25_wp ) )**2
2057            !
2058            zdifml_sc(ji,jj) = pp_dif_ml * phml(ji,jj) * zvel_sc_ml
2059            zvisml_sc(ji,jj) = pp_vis_ml * zdifml_sc(ji,jj)
2060            !
2061            IF ( l_pyc(ji,jj) ) THEN
2062               zdifpyc_n_sc(ji,jj) = pp_dif_pyc * zvel_sc_ml * pdh(ji,jj)
2063               zvispyc_n_sc(ji,jj) = 0.09_wp * zvel_sc_pyc * ( 1.0_wp - phbl(ji,jj) / pdh(ji,jj) )**2 *   &
2064                  &                  ( 0.005_wp  * ( av_u_ml(ji,jj) - av_u_bl(ji,jj) )**2 +     &
2065                  &                    0.0075_wp * ( av_v_ml(ji,jj) - av_v_bl(ji,jj) )**2 ) /   &
2066                  &                  pdh(ji,jj)
2067               zvispyc_n_sc(ji,jj) = pp_vis_pyc * zvel_sc_ml * pdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac
2068               !
2069               IF ( l_shear(ji,jj) .AND. n_ddh(ji,jj) /= 2 ) THEN
2070                  ztmp = pp_vispyc_shr * ( pshear(ji,jj) * phbl(ji,jj) )**pthird * phbl(ji,jj)
2071                  zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + ztmp
2072                  zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + ztmp
2073               ENDIF
2074               !
2075               zdifpyc_s_sc(ji,jj) = pwb_ent(ji,jj) + 0.0025_wp * zvel_sc_pyc * ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) *   &
2076                  &                                   ( av_b_ml(ji,jj) - av_b_bl(ji,jj) )
2077               zvispyc_s_sc(ji,jj) = 0.09_wp * ( pwb_min(ji,jj) + 0.0025_wp * zvel_sc_pyc *                 &
2078                  &                                               ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) *   &
2079                  &                                               ( av_b_ml(ji,jj) - av_b_bl(ji,jj) ) )
2080               zdifpyc_s_sc(ji,jj) = 0.09_wp * zdifpyc_s_sc(ji,jj) * zstab_fac
2081               zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac
2082               !
2083               zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5_wp * zdifpyc_n_sc(ji,jj) )
2084               zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) )
2085               
2086               zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) /   &
2087                  &                           ( zdifml_sc(ji,jj) + epsln ) )**p2third
2088               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 )
2089            ELSE
2090               zdifpyc_n_sc(ji,jj) = pp_dif_pyc * zvel_sc_ml * pdh(ji,jj)   ! ag 19/03
2091               zdifpyc_s_sc(ji,jj) = 0.0_wp   ! ag 19/03
2092               zvispyc_n_sc(ji,jj) = pp_vis_pyc * zvel_sc_ml * pdh(ji,jj)   ! ag 19/03
2093               zvispyc_s_sc(ji,jj) = 0.0_wp   ! ag 19/03
2094               IF(l_coup(ji,jj) ) THEN   ! ag 19/03
2095                  ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub|
2096                  !     already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90
2097                  !  Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub)
2098                  ! wet-cell averaging ..
2099                  zmsku = 0.5_wp * ( 2.0_wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
2100                  zmskv = 0.5_wp * ( 2.0_wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
2101                  zb_coup(ji,jj) = 0.4_wp * SQRT(-1.0_wp * rCdU_bot(ji,jj) *   &
2102                     &             SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2   &
2103                     &                  + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) )
2104                 
2105                  zz_b = -1.0_wp * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)   ! ag 19/03
2106                  zc_coup_vis(ji,jj) = -0.5_wp * ( 0.5_wp * zvisml_sc(ji,jj) / phml(ji,jj) - zb_coup(ji,jj) ) /   &
2107                     &                 ( phml(ji,jj) + zz_b )   ! ag 19/03
2108                  zz_b = -1.0_wp * phml(ji,jj) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)   ! ag 19/03
2109                  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 ) /   &
2110                     &                                  zvisml_sc(ji,jj)   ! ag 19/03
2111                  zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) /   &
2112                     &                           zdifml_sc(ji,jj) )**p2third
2113                  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 +   &
2114                     &                 1.5_wp * ( zdifml_sc(ji,jj) / phml(ji,jj) ) * zbeta_d_sc(ji,jj) *   &
2115                     &                          SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) - zb_coup(ji,jj) ) / zz_b   ! ag 19/03
2116               ELSE   ! ag 19/03
2117                  zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) /   &
2118                     &                           ( zdifml_sc(ji,jj) + epsln ) )**p2third   ! ag 19/03
2119                  zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) /   &
2120                     &                         ( zvisml_sc(ji,jj) + epsln )   ! ag 19/03
2121               ENDIF   ! ag 19/03
2122            ENDIF      ! ag 19/03
2123         ELSE
2124            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)
2125            zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
2126         END IF
2127      END_2D
2128      !
2129      DO_2D( 0, 0, 0, 0 )
2130         IF ( l_conv(ji,jj) ) THEN
2131            DO jk = 2, nmld(ji,jj)   ! Mixed layer diffusivity
2132               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
2133               pdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
2134               pviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_v_sc(ji,jj) * zznd_ml ) *   &
2135                  &                ( 1.0_wp - 0.5_wp * zznd_ml**2 )
2136            END DO
2137            !
2138            ! Coupling to bottom
2139            !
2140            IF ( l_coup(ji,jj) ) THEN                                                         ! ag 19/03
2141               DO jk = mbkt(ji,jj), nmld(ji,jj), -1                                           ! ag 19/03
2142                  zz_b = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )   ! ag 19/03
2143                  pviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2    ! ag 19/03
2144                  pdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2    ! ag 19/03
2145               END DO                                                                         ! ag 19/03
2146            ENDIF                                                                             ! ag 19/03
2147            ! Pycnocline
2148            IF ( l_pyc(ji,jj) ) THEN 
2149               ! Diffusivity and viscosity profiles in the pycnocline given by
2150               ! cubic polynomial. Note, if l_pyc TRUE can't be coupled to seabed.
2151               za_cubic = 0.5_wp
2152               zb_d_cubic = -1.75_wp * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj)
2153               zd_d_cubic = ( pdh(ji,jj) * zdifml_sc(ji,jj) / phml(ji,jj) * SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) *   &
2154                  &           ( 2.5_wp * zbeta_d_sc(ji,jj) - 1.0_wp ) - 0.85_wp * zdifpyc_s_sc(ji,jj) ) /            &
2155                  &           MAX( zdifpyc_n_sc(ji,jj), 1.0e-8_wp )
2156               zd_d_cubic = zd_d_cubic - zb_d_cubic - 2.0_wp * ( 1.0_wp - za_cubic  - zb_d_cubic )
2157               zc_d_cubic = 1.0_wp - za_cubic - zb_d_cubic - zd_d_cubic
2158               zb_v_cubic = -1.75_wp * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)
2159               zd_v_cubic = ( 0.5_wp * zvisml_sc(ji,jj) * pdh(ji,jj) / phml(ji,jj) - 0.85_wp * zvispyc_s_sc(ji,jj) ) /   &
2160                  &           MAX( zvispyc_n_sc(ji,jj), 1.0e-8_wp )
2161               zd_v_cubic = zd_v_cubic - zb_v_cubic - 2.0_wp * ( 1.0_wp - za_cubic - zb_v_cubic )
2162               zc_v_cubic = 1.0_wp - za_cubic - zb_v_cubic - zd_v_cubic
2163               DO jk = nmld(ji,jj) , nbld(ji,jj)
2164                  zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / MAX(pdh(ji,jj), 1.0e-6_wp )
2165                  ztmp = ( 1.75_wp * zznd_pyc - 0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 )
2166                  !
2167                  pdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) *   &
2168                     &                ( za_cubic + zb_d_cubic * zznd_pyc + zc_d_cubic * zznd_pyc**2 + zd_d_cubic * zznd_pyc**3 )
2169                  !
2170                  pdiffut(ji,jj,jk) = pdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ztmp
2171                  pviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) *   &
2172                     &                ( za_cubic + zb_v_cubic * zznd_pyc + zc_v_cubic * zznd_pyc**2 + zd_v_cubic * zznd_pyc**3 )
2173                  pviscos(ji,jj,jk) = pviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ztmp
2174               END DO
2175   !                  IF ( pdhdt(ji,jj) > 0._wp ) THEN
2176   !                     zdiffut(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 )
2177   !                     zviscos(ji,jj,nbld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,nbld(ji,jj)+1,Kmm), 1.0e-6 )
2178   !                  ELSE
2179   !                     zdiffut(ji,jj,nbld(ji,jj)) = 0._wp
2180   !                     zviscos(ji,jj,nbld(ji,jj)) = 0._wp
2181   !                  ENDIF
2182            ENDIF
2183         ELSE
2184            ! Stable conditions
2185            DO jk = 2, nbld(ji,jj)
2186               zznd_ml = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
2187               pdiffut(ji,jj,jk) = 0.75_wp * zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml )**1.5_wp
2188               pviscos(ji,jj,jk) = 0.375_wp * zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml ) * ( 1.0_wp - zznd_ml**2 )
2189            END DO
2190            !
2191            IF ( pdhdt(ji,jj) > 0.0_wp ) THEN
2192               pdiffut(ji,jj,nbld(ji,jj)) = MAX( pdhdt(ji,jj), 1.0e-6_wp) * e3w(ji, jj, nbld(ji,jj), Kmm)
2193               pviscos(ji,jj,nbld(ji,jj)) = pdiffut(ji,jj,nbld(ji,jj))
2194            ENDIF
2195         ENDIF   ! End if ( l_conv )
2196         !
2197      END_2D
2198      IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) )   ! BBL-coupling velocity scale
2199      IF( ln_timing ) CALL timing_stop('zdf_osm_dv')
2200      !
2201   END SUBROUTINE zdf_osm_diffusivity_viscosity
2202
2203   SUBROUTINE zdf_osm_fgr_terms( Kmm, kp_ext, phbl, phml, pdh,                              &
2204      &                          pdhdt, pshear, pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext,   &
2205      &                          pdiffut, pviscos )
2206      !!---------------------------------------------------------------------
2207      !!                 ***  ROUTINE zdf_osm_fgr_terms ***
2208      !!
2209      !! ** Purpose : Compute non-gradient terms in flux-gradient relationship
2210      !!
2211      !! ** Method  :
2212      !!
2213      !!----------------------------------------------------------------------
2214      INTEGER,                          INTENT(in   ) ::   Kmm            ! Time-level index
2215      INTEGER,  DIMENSION(A2D(0)),      INTENT(in   ) ::   kp_ext         ! Offset for external level
2216      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth
2217      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth
2218      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline depth
2219      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency
2220      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pshear         ! Shear production
2221      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients
2222      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients
2223      REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
2224      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(in   ) ::   pdiffut        ! t-diffusivity
2225      REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(in   ) ::   pviscos        ! Viscosity
2226      !
2227      REAL(wp), DIMENSION(A2D(0))     ::   zalpha_pyc   !
2228      REAL(wp), DIMENSION(A2D(0),jpk) ::   zdbdz_pyc    ! Parametrised gradient of buoyancy in the pycnocline
2229      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles
2230      !
2231      INTEGER                     ::   ji, jj, jk, jkm_bld, jkf_mld, jkm_mld   ! Loop indices
2232      INTEGER                     ::   istat                                   ! Memory allocation status
2233      REAL(wp)                    ::   zznd_d, zznd_ml, zznd_pyc, znd          ! Temporary non-dimensional depths
2234      REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_1,zsc_ws_1                      ! Temporary scales
2235      REAL(wp), DIMENSION(A2D(0)) ::   zsc_uw_1, zsc_uw_2                      ! Temporary scales
2236      REAL(wp), DIMENSION(A2D(0)) ::   zsc_vw_1, zsc_vw_2                      ! Temporary scales
2237      REAL(wp), DIMENSION(A2D(0)) ::   ztau_sc_u                               ! Dissipation timescale at base of WML
2238      REAL(wp)                    ::   zbuoy_pyc_sc, zdelta_pyc                !
2239      REAL(wp)                    ::   zl_c,zl_l,zl_eps                        ! Used to calculate turbulence length scale
2240      REAL(wp), DIMENSION(A2D(0)) ::   za_cubic, zb_cubic                      ! Coefficients in cubic polynomial specifying
2241      REAL(wp), DIMENSION(A2D(0)) ::   zc_cubic, zd_cubic                      ! diffusivity in pycnocline
2242      REAL(wp), DIMENSION(A2D(0)) ::   zwt_pyc_sc_1, zws_pyc_sc_1              !
2243      REAL(wp), DIMENSION(A2D(0)) ::   zzeta_pyc                               !
2244      REAL(wp)                    ::   zomega, zvw_max                         !
2245      REAL(wp), DIMENSION(A2D(0)) ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes
2246      REAL(wp), DIMENSION(A2D(0)) ::   zwth_ent,zws_ent                        ! at the top of the pycnocline
2247      REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_pyc, zsc_ws_pyc                 ! Scales for pycnocline transport term
2248      REAL(wp)                    ::   ztmp                                    !
2249      REAL(wp)                    ::   ztgrad, zsgrad, zbgrad                  ! Variables used to calculate pycnocline gradients
2250      REAL(wp)                    ::   zugrad, zvgrad                          ! Variables for calculating pycnocline shear
2251      REAL(wp)                    ::   zdtdz_pyc                               ! Parametrized gradient of temperature in pycnocline
2252      REAL(wp)                    ::   zdsdz_pyc                               ! Parametrised gradient of salinity in pycnocline
2253      REAL(wp)                    ::   zdudz_pyc                               ! u-shear across the pycnocline
2254      REAL(wp)                    ::   zdvdz_pyc                               ! v-shear across the pycnocline
2255      !!----------------------------------------------------------------------
2256      !
2257      IF( ln_timing ) CALL timing_start('zdf_osm_ft')
2258      !
2259      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
2260      !  Pycnocline gradients for scalars and velocity
2261      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2262      CALL zdf_osm_pycnocline_buoyancy_profiles( Kmm, kp_ext, zdbdz_pyc, zalpha_pyc, pdh,    &
2263         &                                       phbl, pdbdz_bl_ext, phml, pdhdt )
2264      !
2265      ! Auxiliary indices
2266      ! -----------------
2267      jkm_bld = 0
2268      jkf_mld = jpk
2269      jkm_mld = 0
2270      DO_2D( 0, 0, 0, 0 )
2271         IF ( nbld(ji,jj) > jkm_bld ) jkm_bld = nbld(ji,jj)
2272         IF ( nmld(ji,jj) < jkf_mld ) jkf_mld = nmld(ji,jj)
2273         IF ( nmld(ji,jj) > jkm_mld ) jkm_mld = nmld(ji,jj)
2274      END_2D
2275      !
2276      ! Stokes term in scalar flux, flux-gradient relationship
2277      ! ------------------------------------------------------
2278      WHERE ( l_conv(A2D(0)) )
2279         zsc_wth_1(:,:) = swstrl(A2D(0))**3 * swth0(A2D(0)) / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )
2280         zsc_ws_1(:,:)  = swstrl(A2D(0))**3 * sws0(A2D(0))  / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )
2281      ELSEWHERE
2282         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D(0))
2283         zsc_ws_1(:,:)  = 2.0_wp * swsav(A2D(0))
2284      ENDWHERE
2285      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
2286         IF ( l_conv(ji,jj) ) THEN
2287            IF ( jk <= nmld(ji,jj) ) THEN
2288               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2289               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) *   &
2290                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj)
2291               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) *   &
2292                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj)
2293            END IF
2294         ELSE   ! Stable conditions
2295            IF ( jk <= nbld(ji,jj) ) THEN
2296               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2297               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) *   &
2298                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj)
2299               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15_wp * EXP( -0.85_wp * zznd_d ) *   &
2300                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj)
2301            END IF
2302         END IF   ! Check on l_conv
2303      END_3D
2304      !
2305      IF ( ln_dia_osm ) THEN
2306         IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
2307         IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
2308      END IF
2309      !
2310      ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use
2311      ! svstr since term needs to go to zero as swstrl goes to zero)
2312      ! ---------------------------------------------------------------------
2313      WHERE ( l_conv(A2D(0)) )
2314         zsc_uw_1(:,:) = ( swstrl(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) /   &
2315            &                                  MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp )
2316         zsc_uw_2(:,:) = ( swstrl(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) /   &
2317            &                                  MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp )
2318         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 ) /   &
2319            &            ( ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln )
2320      ELSEWHERE
2321         zsc_uw_1(:,:) = sustar(A2D(0))**2
2322         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 ) /   &
2323            &            ( svstr(A2D(0))**2 + epsln )
2324      ENDWHERE
2325      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
2326         IF ( l_conv(ji,jj) ) THEN
2327            IF ( jk <= nmld(ji,jj) ) THEN
2328               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2329               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05_wp   * EXP( -0.4_wp * zznd_d ) * zsc_uw_1(ji,jj) +     &
2330                  &                                  0.00125_wp * EXP( -1.0_wp * zznd_d ) * zsc_uw_2(ji,jj) ) *   &
2331                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) )
2332               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65_wp *  0.15_wp * EXP( -1.0_wp * zznd_d ) *                 &
2333                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_vw_1(ji,jj)
2334            END IF
2335         ELSE   ! Stable conditions
2336            IF ( jk <= nbld(ji,jj) ) THEN   ! Corrected to nbld
2337               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2338               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75_wp * 1.3_wp * EXP( -0.5_wp * zznd_d ) *             &
2339                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_uw_1(ji,jj)
2340            END IF
2341         END IF
2342      END_3D
2343      !
2344      ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio
2345      ! (X0.3) and pressure (X0.5)]
2346      ! ----------------------------------------------------------------------
2347      WHERE ( l_conv(A2D(0)) )
2348         zsc_wth_1(:,:) = swbav(A2D(0)) * swth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) /   &
2349            &             ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )
2350         zsc_ws_1(:,:)  = swbav(A2D(0)) * sws0(A2D(0))  * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) /   &
2351            &             ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )
2352      ELSEWHERE
2353         zsc_wth_1(:,:) = 0.0_wp
2354         zsc_ws_1(:,:)  = 0.0_wp
2355      ENDWHERE
2356      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
2357         IF ( l_conv(ji,jj) ) THEN
2358            IF ( jk <= nmld(ji,jj) ) THEN
2359               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
2360               ! Calculate turbulent time scale
2361               zl_c   = 0.9_wp * ( 1.0_wp - EXP( -5.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) *                         &
2362                  &     ( 1.0_wp - EXP( -15.0_wp * ( 1.2_wp - zznd_ml ) ) )
2363               zl_l   = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml + zznd_ml**3 / 3.0_wp ) ) ) *                         &
2364                  &     ( 1.0_wp - EXP( -8.0_wp  * ( 1.15_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) )
2365               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 )
2366               ! Non-gradient buoyancy terms
2367               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 )
2368               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 )
2369            END IF
2370         ELSE   ! Stable conditions
2371            IF ( jk <= nbld(ji,jj) ) THEN
2372               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
2373               ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
2374            END IF
2375         END IF
2376      END_3D
2377      DO_2D( 0, 0, 0, 0 )
2378         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN
2379            ztau_sc_u(ji,jj)    = phml(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *                             &
2380               &                ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )**1.5_wp )
2381            zwth_ent(ji,jj)     = -0.003_wp * ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *   &
2382               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_dt_ml(ji,jj)
2383            zws_ent(ji,jj)      = -0.003_wp * ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *   &
2384               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_ds_ml(ji,jj)
2385            IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) ) THEN
2386               zbuoy_pyc_sc        = 2.0_wp * MAX( av_db_ml(ji,jj), 0.0_wp ) / pdh(ji,jj)
2387               zdelta_pyc          = ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird /   &
2388                  &                       SQRT( MAX( zbuoy_pyc_sc, ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) )
2389               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) ) *   &
2390                  &                     zdelta_pyc**2 / pdh(ji,jj)
2391               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) ) *   &
2392                  &                     zdelta_pyc**2 / pdh(ji,jj)
2393               zzeta_pyc(ji,jj)    = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * shol(ji,jj) ) ) )
2394            END IF
2395         END IF
2396      END_2D
2397      DO_3D( 0, 0, 0, 0, 2, jkm_bld )
2398         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk <= nbld(ji,jj) ) ) THEN
2399            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
2400            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) -                                                                                &
2401               &              0.045_wp * ( ( zwth_ent(ji,jj) * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 &
2402               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp )
2403            ghams(ji,jj,jk) = ghams(ji,jj,jk) -                                                                                &
2404               &              0.045_wp * ( ( zws_ent(ji,jj)  * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 &
2405               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp )
2406            IF ( dh(ji,jj) < 0.2_wp * hbl(ji,jj) .AND. nbld(ji,jj) - nmld(ji,jj) > 3 ) THEN
2407               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp  * zwt_pyc_sc_1(ji,jj) *                              &
2408                  &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        &
2409                  &                                pdh(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird
2410               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp  * zws_pyc_sc_1(ji,jj) *                              &
2411                  &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        &
2412                  &                                pdh(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird
2413            END IF
2414         END IF   ! End of pycnocline
2415      END_3D
2416      !
2417      IF ( ln_dia_osm ) THEN
2418         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! Upward turb. temperature entrainment flux
2419         IF ( iom_use("zws_ent")  ) CALL iom_put( "zws_ent",  tmask(:,:,1)*zws_ent  )   ! Upward turb. salinity    entrainment flux
2420      END IF
2421      !
2422      zsc_vw_1(:,:) = 0.0_wp
2423      WHERE ( l_conv(A2D(0)) )
2424         zsc_uw_1(:,:) = -1.0_wp * swb0(A2D(0)) * sustar(A2D(0))**2 * phml(A2D(0)) /   &
2425            &            ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )
2426         zsc_uw_2(:,:) =           swb0(A2D(0)) * sustke(A2D(0))    * phml(A2D(0)) /   &
2427            &            ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp )
2428      ELSEWHERE
2429         zsc_uw_1(:,:) = 0.0_wp
2430      ENDWHERE
2431      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
2432         IF ( l_conv(ji,jj) ) THEN
2433            IF ( jk <= nmld(ji,jj) ) THEN
2434               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2435               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3_wp * 0.5_wp *   &
2436                  &                                ( zsc_uw_1(ji,jj) + 0.125_wp * EXP( -0.5_wp * zznd_d ) *       &
2437                  &                                  (   1.0_wp - EXP( -0.5_wp * zznd_d ) ) * zsc_uw_2(ji,jj) )
2438               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
2439            END IF
2440         ELSE   ! Stable conditions
2441            IF ( jk <= nbld(ji,jj) ) THEN
2442               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
2443               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
2444            END IF
2445         ENDIF
2446      END_3D
2447      !
2448      DO_2D( 0, 0, 0, 0 )
2449         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN
2450            IF ( n_ddh(ji,jj) == 0 ) THEN
2451               ! Place holding code. Parametrization needs checking for these conditions.
2452               zomega = ( 0.15_wp * swstrl(ji,jj)**3 + swstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird
2453               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_du_ml(ji,jj)
2454               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_dv_ml(ji,jj)
2455            ELSE
2456               zomega = ( 0.15_wp * swstrl(ji,jj)**3 + swstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird
2457               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_du_ml(ji,jj)
2458               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_dv_ml(ji,jj)
2459            ENDIF
2460            zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * suw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj)
2461            za_cubic(ji,jj) = zuw_bse(ji,jj) - zb_cubic(ji,jj)
2462            zvw_max = 0.7_wp * ff_t(ji,jj) * ( sustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * sustar(ji,jj) * phml(ji,jj) )
2463            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)
2464            zc_cubic(ji,jj) = zvw_bse(ji,jj) - zd_cubic(ji,jj)
2465         END IF
2466      END_2D
2467      DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array.
2468         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN
2469            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
2470            ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse(ji,jj) *                 &
2471               &                                ( za_cubic(ji,jj) * zznd_pyc**2 + zb_cubic(ji,jj) * zznd_pyc**3 ) *   &
2472               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
2473            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse(ji,jj) *                 &
2474               &                                ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) *   &
2475               &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
2476         END IF   ! l_conv .AND. l_pyc
2477      END_3D
2478      !
2479      IF ( ln_dia_osm ) THEN
2480         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0",    wmask*ghamu           )
2481         IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
2482      END IF
2483      !
2484      ! Transport term in flux-gradient relationship [note : includes ROI ratio
2485      ! (X0.3) ]
2486      ! -----------------------------------------------------------------------
2487      WHERE ( l_conv(A2D(0)) )
2488         zsc_wth_1(:,:) = swth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) )
2489         zsc_ws_1(:,:)  = sws0(A2D(0))  / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) )
2490         WHERE ( l_pyc(A2D(0)) )   ! Pycnocline scales
2491            zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_dt_ml(A2D(0))
2492            zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_ds_ml(A2D(0))
2493         END WHERE
2494      ELSEWHERE
2495         zsc_wth_1(:,:) = 2.0 * swthav(A2D(0))
2496         zsc_ws_1(:,:)  =       sws0(A2D(0))
2497      END WHERE
2498      DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) )
2499         IF ( l_conv(ji,jj) ) THEN
2500            IF ( ( jk > 1 ) .AND. ( jk <= nmld(ji,jj) ) ) THEN
2501               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
2502               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 ) ) ) *   &
2503                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
2504               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 ) ) ) *   &
2505                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
2506            END IF
2507            !
2508            ! may need to comment out lpyc block
2509            IF ( l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN   ! Pycnocline
2510               zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
2511               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 ) )
2512               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 ) )
2513            END IF
2514         ELSE
2515            IF( pdhdt(ji,jj) > 0. ) THEN
2516               IF ( ( jk > 1 ) .AND. ( jk <= nbld(ji,jj) ) ) THEN
2517                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2518                  znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
2519                  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 ) ) +   &
2520                     7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_wth_1(ji,jj)
2521                  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 ) ) +   &
2522                     7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_ws_1(ji,jj)
2523               END IF
2524            ENDIF
2525         ENDIF
2526      END_3D
2527      !
2528      WHERE ( l_conv(A2D(0)) )
2529         zsc_uw_1(:,:) = sustar(A2D(0))**2
2530         zsc_vw_1(:,:) = ff_t(A2D(0)) * sustke(A2D(0)) * phml(A2D(0))
2531      ELSEWHERE
2532         zsc_uw_1(:,:) = sustar(A2D(0))**2
2533         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 ) ) *   &
2534            &            zsc_uw_1(:,:)
2535         zsc_vw_1(:,:) = ff_t(A2D(0)) * sustke(A2D(0)) * phbl(A2D(0))
2536         zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) *   &
2537            &            zsc_vw_1(:,:)
2538      ENDWHERE
2539      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
2540         IF ( l_conv(ji,jj) ) THEN
2541            IF ( jk <= nmld(ji,jj) ) THEN
2542               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
2543               zznd_d  = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2544               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +   &
2545                  &              0.3_wp * ( -2.0_wp + 2.5_wp * ( 1.0_wp + 0.1_wp * zznd_ml**4 ) - EXP( -8.0_wp * zznd_ml ) ) *   &
2546                  &              zsc_uw_1(ji,jj)
2547               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) +   &
2548                  &              0.3_wp * 0.1_wp * ( EXP( -1.0_wp * zznd_d ) + EXP( -5.0_wp * ( 1.0_wp - zznd_ml ) ) ) *   &
2549                  &              zsc_vw_1(ji,jj)
2550            END IF
2551         ELSE
2552            IF ( jk <= nbld(ji,jj) ) THEN
2553               znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
2554               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
2555               IF ( zznd_d <= 2.0 ) THEN
2556                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp *                                              &
2557                     &                                ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * zznd_d ) ) *   &
2558                     &                                  ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) ) * zsc_uw_1(ji,jj)
2559               ELSE
2560                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp *   &
2561                     &                                ( 1.0_wp - EXP( -5.0_wp * ( 1.0_wp - znd ) ) ) * zsc_uw_2(ji,jj)
2562               ENDIF
2563               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * SIN( 3.14159_wp * ( 0.65_wp * zznd_d ) ) *   &
2564                  &                                EXP( -0.25_wp * zznd_d**2 ) * zsc_vw_1(ji,jj)
2565               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
2566            END IF
2567         END IF
2568      END_3D
2569      !
2570      IF ( ln_dia_osm ) THEN
2571         IF ( iom_use("ghamu_f") )    CALL iom_put( "ghamu_f",    wmask       *ghamu    )
2572         IF ( iom_use("ghamv_f") )    CALL iom_put( "ghamv_f",    wmask       *ghamv    )
2573         IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
2574         IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
2575         IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
2576         IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
2577      END IF
2578      !
2579      ! Make surface forced velocity non-gradient terms go to zero at the base
2580      ! of the mixed layer.
2581      !
2582      ! Make surface forced velocity non-gradient terms go to zero at the base
2583      ! of the boundary layer.
2584      DO_3D( 0, 0, 0, 0, 2, jkm_bld )
2585         IF ( ( .NOT. l_conv(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN
2586            znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj)   ! ALMG to think about
2587            IF ( znd >= 0.0_wp ) THEN
2588               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) )
2589               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) )
2590            ELSE
2591               ghamu(ji,jj,jk) = 0.0_wp
2592               ghamv(ji,jj,jk) = 0.0_wp
2593            ENDIF
2594         END IF
2595      END_3D
2596      !
2597      ! Pynocline contributions
2598      !
2599      IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Allocate arrays for output of pycnocline gradient/shear profiles
2600         ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat )
2601         IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' )
2602         z3ddz_pyc_1(:,:,:) = 0.0_wp
2603         z3ddz_pyc_2(:,:,:) = 0.0_wp
2604      END IF
2605      DO_3D( 0, 0, 0, 0, 2, jkm_bld )
2606         IF ( l_conv (ji,jj) ) THEN
2607            ! Unstable conditions. Shouldn;t be needed with no pycnocline code.
2608            !                  zugrad = 0.7 * av_du_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
2609            !                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
2610            !                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
2611            !Alan is this right?
2612            !                  zvgrad = ( 0.7 * av_dv_ml(ji,jj) + &
2613            !                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
2614            !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
2615            !                       &      )/ (zdh(ji,jj)  + epsln )
2616            !                  DO jk = 2, nbld(ji,jj) - 1 + ibld_ext
2617            !                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
2618            !                     IF ( znd <= 0.0 ) THEN
2619            !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
2620            !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
2621            !                     ELSE
2622            !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
2623            !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
2624            !                     ENDIF
2625            !                  END DO
2626         ELSE   ! Stable conditions
2627            IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN
2628               ! Pycnocline profile only defined when depth steady of increasing.
2629               IF ( pdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady.
2630                  IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN
2631                     IF ( shol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline
2632                        ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln )
2633                        ztgrad = av_dt_bl(ji,jj) * ztmp
2634                        zsgrad = av_ds_bl(ji,jj) * ztmp
2635                        zbgrad = av_db_bl(ji,jj) * ztmp
2636                        IF ( jk <= nbld(ji,jj) ) THEN
2637                           znd = gdepw(ji,jj,jk,Kmm) * ztmp
2638                           zdtdz_pyc =  ztgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 )
2639                           zdsdz_pyc =  zsgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 )
2640                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc
2641                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc
2642                           IF ( ln_dia_pyc_scl ) THEN
2643                              z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc
2644                              z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc
2645                           END IF
2646                        END IF
2647                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
2648                        ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln )
2649                        ztgrad = av_dt_bl(ji,jj) * ztmp
2650                        zsgrad = av_ds_bl(ji,jj) * ztmp
2651                        zbgrad = av_db_bl(ji,jj) * ztmp
2652                        IF ( jk <= nbld(ji,jj) ) THEN
2653                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp
2654                           zdtdz_pyc =  ztgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 )
2655                           zdsdz_pyc =  zsgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 )
2656                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc
2657                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(ji,jj,jk) * zdsdz_pyc
2658                           IF ( ln_dia_pyc_scl ) THEN
2659                              z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc
2660                              z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc
2661                           END IF
2662                        END IF
2663                     ENDIF   ! IF (shol >=0.5)
2664                  ENDIF      ! IF (av_db_bl> 0.)
2665               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
2666            END IF
2667         END IF
2668      END_3D
2669      IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles
2670         IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )
2671         IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )
2672      END IF
2673      DO_3D( 0, 0, 0, 0, 2, jkm_bld )
2674         IF ( .NOT. l_conv (ji,jj) ) THEN
2675            IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN
2676               zugrad = 3.25_wp * av_du_bl(ji,jj) / phbl(ji,jj)
2677               zvgrad = 2.75_wp * av_dv_bl(ji,jj) / phbl(ji,jj)
2678               IF ( jk <= nbld(ji,jj) ) THEN
2679                  znd = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
2680                  IF ( znd < 1.0 ) THEN
2681                     zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 )
2682                  ELSE
2683                     zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 )
2684                  ENDIF
2685                  zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 )
2686                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + pviscos(ji,jj,jk) * zdudz_pyc
2687                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + pviscos(ji,jj,jk) * zdvdz_pyc
2688                  IF ( ln_dia_pyc_shr ) THEN
2689                     z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc
2690                     z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc
2691                  END IF
2692               END IF
2693            END IF
2694         END IF
2695      END_3D
2696      IF ( ln_dia_pyc_shr ) THEN   ! Output of pycnocline shear profiles
2697         IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )
2698         IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )
2699      END IF
2700      IF ( ln_dia_osm ) THEN
2701         IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
2702         IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
2703      END IF
2704      IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Deallocate arrays used for output of pycnocline gradient/shear profiles
2705         DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 )
2706      END IF
2707      !
2708      DO_2D( 0, 0, 0, 0 )
2709         ghamt(ji,jj,nbld(ji,jj)) = 0.0_wp
2710         ghams(ji,jj,nbld(ji,jj)) = 0.0_wp
2711         ghamu(ji,jj,nbld(ji,jj)) = 0.0_wp
2712         ghamv(ji,jj,nbld(ji,jj)) = 0.0_wp
2713      END_2D
2714      !
2715      IF ( ln_dia_osm ) THEN
2716         IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
2717         IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
2718         IF ( iom_use("zviscos") ) THEN
2719            osmdia3d(A2D(0),:) = wmask(A2D(0),:) * pviscos; CALL iom_put( "zviscos", osmdia3d )
2720         END IF
2721      END IF
2722      !
2723      IF( ln_timing ) CALL timing_stop('zdf_osm_ft')
2724      !
2725   END SUBROUTINE zdf_osm_fgr_terms
2726
2727   SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm, pmld, pdtdx, pdtdy, pdsdx,   &
2728      &                                          pdsdy, pdbdx_mle, pdbdy_mle, pdbds_mle )
2729      !!----------------------------------------------------------------------
2730      !!          ***  ROUTINE zdf_osm_zmld_horizontal_gradients  ***
2731      !!
2732      !! ** Purpose : Calculates horizontal gradients of buoyancy for use with
2733      !!              Fox-Kemper parametrization
2734      !!
2735      !! ** Method  :
2736      !!
2737      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
2738      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
2739      !!
2740      !!----------------------------------------------------------------------
2741      INTEGER,                      INTENT(in   ) ::   Kmm          ! Time-level index
2742      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pmld         ! == Estimated FK BLD used for MLE horizontal gradients == !
2743      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdtdx        ! Horizontal gradient for Fox-Kemper parametrization
2744      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdtdy        ! Horizontal gradient for Fox-Kemper parametrization
2745      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdsdx        ! Horizontal gradient for Fox-Kemper parametrization
2746      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdsdy        ! Horizontal gradient for Fox-Kemper parametrization
2747      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdx_mle    ! MLE horiz gradients at u points
2748      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdy_mle    ! MLE horiz gradients at v points
2749      REAL(wp), DIMENSION(A2D(0)),  INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient
2750      !
2751      ! Local variables
2752      INTEGER                          ::   ji, jj, jk   ! Dummy loop indices
2753      INTEGER                          ::   ikt, ikmax   ! Local integers     
2754      REAL(wp)                         ::   zc
2755      REAL(wp)                         ::   zN2_c        ! Local buoyancy difference from 10m value
2756      REAL(wp), DIMENSION(A2D(1))      ::   ztm
2757      REAL(wp), DIMENSION(A2D(1))      ::   zsm
2758      REAL(wp), DIMENSION(A2D(1),jpts) ::   ztsm_midu
2759      REAL(wp), DIMENSION(A2D(1),jpts) ::   ztsm_midv
2760      REAL(wp), DIMENSION(A2D(1),jpts) ::   zabu
2761      REAL(wp), DIMENSION(A2D(1),jpts) ::   zabv
2762      REAL(wp), DIMENSION(A2D(1))      ::   zmld_midu
2763      REAL(wp), DIMENSION(A2D(1))      ::   zmld_midv
2764      !
2765      IF( ln_timing ) CALL timing_start('zdf_osm_zhg')
2766      !
2767      ! ==  MLD used for MLE  ==!
2768      mld_prof(:,:) = nlb10   ! Initialization to the number of w ocean point
2769      pmld(:,:)  = 0.0_wp     ! Here hmlp used as a dummy variable, integrating vertically N^2
2770      zN2_c = grav * rn_osm_mle_rho_c * r1_rho0   ! Convert density criteria into N^2 criteria
2771      DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )
2772         ikt = mbkt(ji,jj)
2773         pmld(ji,jj) = pmld(ji,jj) + MAX( rn2b(ji,jj,jk), 0.0_wp ) * e3w(ji,jj,jk,Kmm)
2774         IF( pmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
2775      END_3D
2776      DO_2D( 1, 1, 1, 1 )
2777         mld_prof(ji,jj) = MAX( mld_prof(ji,jj), nbld(ji,jj) )   ! Ensure mld_prof .ge. nbld
2778         pmld(ji,jj)     = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
2779      END_2D
2780      !
2781      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )   ! Max level of the computation
2782      ztm(:,:) = 0.0_wp
2783      zsm(:,:) = 0.0_wp
2784      DO_3D( 1, 1, 1, 1, 1, ikmax )
2785         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj) - jk ), 1  ), KIND=wp )   ! zc being 0 outside the ML t-points
2786         ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm)
2787         zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm)
2788      END_3D
2789      ! Average temperature and salinity
2790      ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), pmld(:,:) )
2791      zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), pmld(:,:) )
2792      ! Calculate horizontal gradients at u & v points
2793      zmld_midu(:,:)   =  0.0_wp
2794      ztsm_midu(:,:,:) = 10.0_wp
2795      DO_2D( 1, 0, 0, 0 )
2796         pdtdx(ji,jj)            = ( ztm(ji+1,jj) - ztm(ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
2797         pdsdx(ji,jj)            = ( zsm(ji+1,jj) - zsm(ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
2798         zmld_midu(ji,jj)        = 0.25_wp * ( pmld(ji+1,jj) + pmld(ji,jj))
2799         ztsm_midu(ji,jj,jp_tem) =  0.5_wp * ( ztm( ji+1,jj)  + ztm( ji,jj) )
2800         ztsm_midu(ji,jj,jp_sal) =  0.5_wp * ( zsm( ji+1,jj)  + zsm( ji,jj) )
2801      END_2D
2802      zmld_midv(:,:)   =  0.0_wp
2803      ztsm_midv(:,:,:) = 10.0_wp
2804      DO_2D( 0, 0, 1, 0 )
2805         pdtdy(ji,jj)            = ( ztm(ji,jj+1) - ztm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
2806         pdsdy(ji,jj)            = ( zsm(ji,jj+1) - zsm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
2807         zmld_midv(ji,jj)        = 0.25_wp * ( pmld(ji,jj+1) + pmld( ji,jj) )
2808         ztsm_midv(ji,jj,jp_tem) =  0.5_wp * ( ztm( ji,jj+1)  + ztm( ji,jj) )
2809         ztsm_midv(ji,jj,jp_sal) =  0.5_wp * ( zsm( ji,jj+1)  + zsm( ji,jj) )
2810      END_2D
2811      CALL eos_rab( ztsm_midu, zmld_midu, zabu, Kmm )
2812      CALL eos_rab( ztsm_midv, zmld_midv, zabv, Kmm )
2813      DO_2D( 1, 0, 0, 0 )
2814         pdbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) )
2815      END_2D
2816      DO_2D( 0, 0, 1, 0 )
2817         pdbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) )
2818      END_2D
2819      DO_2D( 0, 0, 0, 0 )
2820         pdbds_mle(ji,jj) = SQRT( 0.5_wp * ( pdbdx_mle(ji,  jj) * pdbdx_mle(ji,  jj) + pdbdy_mle(ji,jj  ) * pdbdy_mle(ji,jj  ) +   &
2821            &                                pdbdx_mle(ji-1,jj) * pdbdx_mle(ji-1,jj) + pdbdy_mle(ji,jj-1) * pdbdy_mle(ji,jj-1) ) )
2822      END_2D
2823      !
2824      IF( ln_timing ) CALL timing_stop('zdf_osm_zhg')
2825      !
2826   END SUBROUTINE zdf_osm_zmld_horizontal_gradients
2827
2828   SUBROUTINE zdf_osm_osbl_state_fk( Kmm, pwb_fk, phbl, phmle, pwb_ent,   &
2829      &                              pdbds_mle )
2830      !!---------------------------------------------------------------------
2831      !!               ***  ROUTINE zdf_osm_osbl_state_fk  ***
2832      !!
2833      !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is
2834      !!              returned in the logicals l_pyc, l_flux and ldmle. Used
2835      !!              with Fox-Kemper scheme.
2836      !!                l_pyc  :: determines whether pycnocline flux-grad
2837      !!                          relationship needs to be determined
2838      !!                l_flux :: determines whether effects of surface flux
2839      !!                          extend below the base of the OSBL
2840      !!                ldmle  :: determines whether the layer with MLE is
2841      !!                          increasing with time or if base is relaxing
2842      !!                          towards hbl
2843      !!
2844      !! ** Method  :
2845      !!
2846      !!----------------------------------------------------------------------     
2847      ! Outputs
2848      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index
2849      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pwb_fk
2850      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl        ! BL depth
2851      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phmle       ! MLE depth
2852      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux
2853      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient
2854      !
2855      ! Local variables
2856      INTEGER                     ::   ji, jj, jk        ! Dummy loop indices
2857      REAL(wp), DIMENSION(A2D(0)) ::   znd_param
2858      REAL(wp)                    ::   zthermal, zbeta
2859      REAL(wp)                    ::   zbuoy
2860      REAL(wp)                    ::   ztmp
2861      REAL(wp)                    ::   zpe_mle_layer
2862      REAL(wp)                    ::   zpe_mle_ref
2863      REAL(wp)                    ::   zdbdz_mle_int
2864      !
2865      IF( ln_timing ) CALL timing_start('zdf_osm_osf')
2866      !
2867      znd_param(A2D(0)) = 0.0_wp
2868      !
2869      DO_2D( 0, 0, 0, 0 )
2870         ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
2871         pwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * pdbds_mle(ji,jj) * pdbds_mle(ji,jj)
2872      END_2D
2873      !
2874      DO_2D( 0, 0, 0, 0 )
2875         !
2876         IF ( l_conv(ji,jj) ) THEN
2877            IF ( phmle(ji,jj) > 1.2_wp * phbl(ji,jj) ) THEN
2878               av_t_mle(ji,jj) = ( av_t_mle(ji,jj) * phmle(ji,jj) - av_t_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) )
2879               av_s_mle(ji,jj) = ( av_s_mle(ji,jj) * phmle(ji,jj) - av_s_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) )
2880               av_b_mle(ji,jj) = ( av_b_mle(ji,jj) * phmle(ji,jj) - av_b_bl(ji,jj) * phbl(ji,jj) ) / ( phmle(ji,jj) - phbl(ji,jj) )
2881               zdbdz_mle_int = ( av_b_bl(ji,jj) - ( 2.0_wp * av_b_mle(ji,jj) - av_b_bl(ji,jj) ) ) / ( phmle(ji,jj) - phbl(ji,jj) )
2882               ! Calculate potential energies of actual profile and reference profile
2883               zpe_mle_layer = 0.0_wp
2884               zpe_mle_ref   = 0.0_wp
2885               zthermal = rab_n(ji,jj,1,jp_tem)
2886               zbeta    = rab_n(ji,jj,1,jp_sal)
2887               DO jk = nbld(ji,jj), mld_prof(ji,jj)
2888                  zbuoy         = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) )
2889                  zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
2890                  zpe_mle_ref   = zpe_mle_ref   + ( av_b_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) ) *   &
2891                     &                            gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
2892               END DO
2893               ! Non-dimensional parameter to diagnose the presence of thermocline
2894               znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) /   &
2895                  &               ( MAX( pwb_fk(ji,jj), 1e-10 ) * phmle(ji,jj) )
2896            END IF
2897         END IF
2898         !
2899      END_2D
2900      !
2901      ! Diagnosis
2902      DO_2D( 0, 0, 0, 0 )
2903         !
2904         IF ( l_conv(ji,jj) ) THEN
2905            IF ( -2.0_wp * pwb_fk(ji,jj) / pwb_ent(ji,jj) > 0.5_wp ) THEN
2906               IF ( phmle(ji,jj) > 1.2_wp * phbl(ji,jj) ) THEN   ! MLE layer growing
2907                  IF ( znd_param (ji,jj) > 100.0_wp ) THEN   ! Thermocline present
2908                     l_flux(ji,jj) = .FALSE.
2909                     l_mle(ji,jj)  = .FALSE.
2910                  ELSE   ! Thermocline not present
2911                     l_flux(ji,jj) = .TRUE.
2912                     l_mle(ji,jj)  = .TRUE.
2913                  ENDIF  ! znd_param > 100
2914                  !
2915                  IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh ) THEN
2916                     l_pyc(ji,jj) = .FALSE.
2917                  ELSE
2918                     l_pyc(ji,jj) = .TRUE.
2919                  ENDIF
2920               ELSE   ! MLE layer restricted to OSBL or just below
2921                  IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh ) THEN   ! Weak stratification MLE layer can grow
2922                     l_pyc(ji,jj)  = .FALSE.
2923                     l_flux(ji,jj) = .TRUE.
2924                     l_mle(ji,jj)  = .TRUE.
2925                  ELSE   ! Strong stratification
2926                     l_pyc(ji,jj)  = .TRUE.
2927                     l_flux(ji,jj) = .FALSE.
2928                     l_mle(ji,jj)  = .FALSE.
2929                  END IF   ! av_db_bl < rn_mle_thresh_bl and
2930               END IF   ! phmle > 1.2 phbl
2931            ELSE
2932               l_pyc(ji,jj)  = .TRUE.
2933               l_flux(ji,jj) = .FALSE.
2934               l_mle(ji,jj)  = .FALSE.
2935               IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE.
2936            END IF   !  -2.0 * pwb_fk(ji,jj) / pwb_ent > 0.5
2937         ELSE   ! Stable Boundary Layer
2938            l_pyc(ji,jj)  = .FALSE.
2939            l_flux(ji,jj) = .FALSE.
2940            l_mle(ji,jj)  = .FALSE.
2941         END IF   ! l_conv
2942         !
2943      END_2D
2944      !
2945      IF( ln_timing ) CALL timing_stop('zdf_osm_osf')
2946      !
2947   END SUBROUTINE zdf_osm_osbl_state_fk
2948
2949   SUBROUTINE zdf_osm_mle_parameters( Kmm, kmld_prof, pmld, phmle, pvel_mle,   &
2950      &                               pdiff_mle, pdbds_mle, phbl, pwb0tot )
2951      !!----------------------------------------------------------------------
2952      !!               ***  ROUTINE zdf_osm_mle_parameters  ***
2953      !!
2954      !! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates
2955      !!              the mixed layer eddy fluxes for buoyancy, heat and
2956      !!              salinity.
2957      !!
2958      !! ** Method  :
2959      !!
2960      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
2961      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
2962      !!
2963      !!----------------------------------------------------------------------
2964      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index
2965      INTEGER,  DIMENSION(:,:),    INTENT(inout) ::   kmld_prof
2966      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == !
2967      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phmle       ! MLE depth
2968      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pvel_mle    ! Velocity scale for dhdt with stable ML and FK
2969      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdiff_mle   ! Extra MLE vertical diff
2970      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient
2971      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl        ! BL depth
2972      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb0tot     ! Total surface buoyancy flux including insolation
2973      !
2974      ! Local variables
2975      INTEGER  ::   ji, jj, jk   ! Dummy loop indices
2976      REAL(wp) ::   ztmp
2977      REAL(wp) ::   zdbdz
2978      REAL(wp) ::   zdtdz
2979      REAL(wp) ::   zdsdz
2980      REAL(wp) ::   zthermal
2981      REAL(wp) ::   zbeta
2982      REAL(wp) ::   zbuoy
2983      REAL(wp) ::   zdb_mle
2984      !
2985      IF( ln_timing ) CALL timing_start('zdf_osm_mp')
2986      !
2987      ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE
2988      DO_2D( 0, 0, 0, 0 )
2989         IF ( l_conv(ji,jj) ) THEN
2990            ztmp =  r1_ft(ji,jj) * MIN( 111e3_wp, e1u(ji,jj) ) / rn_osm_mle_lf
2991            ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt
2992            pvel_mle(ji,jj)  = pdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1)
2993            pdiff_mle(ji,jj) = 5e-4_wp * rn_osm_mle_ce * ztmp * pdbds_mle(ji,jj) * phmle(ji,jj)**2
2994         END IF
2995      END_2D
2996      ! Timestep mixed layer eddy depth
2997      DO_2D( 0, 0, 0, 0 )
2998         IF ( l_mle(ji,jj) ) THEN   ! MLE layer growing
2999            ! Buoyancy gradient at base of MLE layer
3000            zthermal = rab_n(ji,jj,1,jp_tem)
3001            zbeta    = rab_n(ji,jj,1,jp_sal)
3002            zbuoy = grav * ( zthermal * ts(ji,jj,kmld_prof(ji,jj)+2,jp_tem,Kmm) -   &
3003               &             zbeta    * ts(ji,jj,kmld_prof(ji,jj)+2,jp_sal,Kmm) )
3004            zdb_mle = av_b_bl(ji,jj) - zbuoy
3005            ! Timestep hmle
3006            hmle(ji,jj) = hmle(ji,jj) + pwb0tot(ji,jj) * rn_Dt / zdb_mle
3007         ELSE
3008            IF ( phmle(ji,jj) > phbl(ji,jj) ) THEN
3009               hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau
3010            ELSE
3011               hmle(ji,jj) = hmle(ji,jj) - 10.0_wp * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau
3012            END IF
3013         END IF
3014         hmle(ji,jj) = MAX( MIN( hmle(ji,jj), ht(ji,jj) ), gdepw(ji,jj,4,Kmm) )
3015         IF ( ln_osm_hmle_limit ) hmle(ji,jj) = MIN( hmle(ji,jj), rn_osm_hmle_limit*hbl(ji,jj) )
3016         hmle(ji,jj) = pmld(ji,jj)   ! For now try just set hmle to pmld
3017      END_2D
3018      !
3019      kmld_prof(:,:) = 4
3020      DO_3D( 0, 0, 0, 0, 5, jpkm1 )
3021         IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) kmld_prof(ji,jj) = MIN( mbkt(ji,jj), jk )
3022      END_3D
3023      DO_2D( 0, 0, 0, 0 )
3024         phmle(ji,jj) = gdepw(ji,jj,kmld_prof(ji,jj),Kmm)
3025      END_2D
3026      !
3027      IF( ln_timing ) CALL timing_stop('zdf_osm_mp')
3028      !
3029   END SUBROUTINE zdf_osm_mle_parameters
3030
3031   SUBROUTINE zdf_osm_init( Kmm )
3032      !!----------------------------------------------------------------------
3033      !!                  ***  ROUTINE zdf_osm_init  ***
3034      !!
3035      !! ** Purpose :   Initialization of the vertical eddy diffivity and
3036      !!      viscosity when using a osm turbulent closure scheme
3037      !!
3038      !! ** Method  :   Read the namosm namelist and check the parameters
3039      !!      called at the first timestep (nit000)
3040      !!
3041      !! ** input   :   Namlists namzdf_osm and namosm_mle
3042      !!
3043      !!----------------------------------------------------------------------
3044      INTEGER, INTENT(in   ) ::   Kmm   ! Time level
3045      !
3046      ! Local variables
3047      INTEGER  ::   ios            ! Local integer
3048      INTEGER  ::   ji, jj, jk     ! Dummy loop indices
3049      REAL(wp) ::   z1_t2
3050      !
3051      REAL(wp), PARAMETER ::   pp_large = -1e10_wp
3052      !
3053      NAMELIST/namzdf_osm/ ln_use_osm_la,    rn_osm_la,      rn_osm_dstokes,      nn_ave,                nn_osm_wave,        &
3054         &                 ln_dia_osm,       rn_osm_hbl0,    rn_zdfosm_adjust_sd, ln_kpprimix,           rn_riinfty,         &
3055         &                 rn_difri,         ln_convmix,     rn_difconv,          nn_osm_wave,           nn_osm_SD_reduce,   &
3056         &                 ln_osm_mle,       rn_osm_hblfrac, rn_osm_bl_thresh,    ln_zdfosm_ice_shelter
3057      ! Namelist for Fox-Kemper parametrization
3058      NAMELIST/namosm_mle/ nn_osm_mle,       rn_osm_mle_ce,     rn_osm_mle_lf,  rn_osm_mle_time,  rn_osm_mle_lat,   &
3059         &                 rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit
3060      !
3061      IF( ln_timing ) CALL timing_start('zdf_osm_init')
3062      READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
3063901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
3064
3065      READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
3066902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
3067      IF(lwm) WRITE ( numond, namzdf_osm )
3068
3069      IF(lwp) THEN                    ! Control print
3070         WRITE(numout,*)
3071         WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
3072         WRITE(numout,*) '~~~~~~~~~~~~'
3073         WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters'
3074         WRITE(numout,*) '      Use  rn_osm_la                                     ln_use_osm_la         = ', ln_use_osm_la
3075         WRITE(numout,*) '      Use  MLE in OBL, i.e. Fox-Kemper param             ln_osm_mle            = ', ln_osm_mle
3076         WRITE(numout,*) '      Turbulent Langmuir number                          rn_osm_la             = ', rn_osm_la
3077         WRITE(numout,*) '      Stokes drift reduction factor                      rn_zdfosm_adjust_sd   = ', rn_zdfosm_adjust_sd
3078         WRITE(numout,*) '      Initial hbl for 1D runs                            rn_osm_hbl0           = ', rn_osm_hbl0
3079         WRITE(numout,*) '      Depth scale of Stokes drift                        rn_osm_dstokes        = ', rn_osm_dstokes
3080         WRITE(numout,*) '      Horizontal average flag                            nn_ave                = ', nn_ave
3081         WRITE(numout,*) '      Stokes drift                                       nn_osm_wave           = ', nn_osm_wave
3082         SELECT CASE (nn_osm_wave)
3083         CASE(0)
3084            WRITE(numout,*) '      Calculated assuming constant La#=0.3'
3085         CASE(1)
3086            WRITE(numout,*) '      Calculated from Pierson Moskowitz wind-waves'
3087         CASE(2)
3088            WRITE(numout,*) '      Calculated from ECMWF wave fields'
3089         END SELECT
3090         WRITE(numout,*) '      Stokes drift reduction                             nn_osm_SD_reduce      = ', nn_osm_SD_reduce
3091         WRITE(numout,*) '      Fraction of hbl to average SD over/fit'
3092         WRITE(numout,*) '      Exponential with nn_osm_SD_reduce = 1 or 2         rn_osm_hblfrac        = ', rn_osm_hblfrac
3093         SELECT CASE (nn_osm_SD_reduce)
3094         CASE(0)
3095            WRITE(numout,*) '     No reduction'
3096         CASE(1)
3097            WRITE(numout,*) '     Average SD over upper rn_osm_hblfrac of BL'
3098         CASE(2)
3099            WRITE(numout,*) '     Fit exponential to slope rn_osm_hblfrac of BL'
3100         END SELECT
3101         WRITE(numout,*) '     Reduce surface SD and depth scale under ice         ln_zdfosm_ice_shelter = ', ln_zdfosm_ice_shelter
3102         WRITE(numout,*) '     Output osm diagnostics                              ln_dia_osm            = ', ln_dia_osm
3103         WRITE(numout,*) '         Threshold used to define BL                     rn_osm_bl_thresh      = ', rn_osm_bl_thresh,   &
3104            &            'm^2/s'
3105         WRITE(numout,*) '     Use KPP-style shear instability mixing              ln_kpprimix           = ', ln_kpprimix
3106         WRITE(numout,*) '     Local Richardson Number limit for shear instability rn_riinfty            = ', rn_riinfty
3107         WRITE(numout,*) '     Maximum shear diffusivity at Rig = 0 (m2/s)         rn_difri              = ', rn_difri
3108         WRITE(numout,*) '     Use large mixing below BL when unstable             ln_convmix            = ', ln_convmix
3109         WRITE(numout,*) '     Diffusivity when unstable below BL (m2/s)           rn_difconv            = ', rn_difconv
3110      ENDIF
3111      !
3112      !                              ! Check wave coupling settings !
3113      !                         ! Further work needed - see ticket #2447 !
3114      IF ( nn_osm_wave == 2 ) THEN
3115         IF (.NOT. ( ln_wave .AND. ln_sdw )) &
3116            & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' )
3117      END IF
3118      !
3119      ! Flags associated with diagnostic output
3120      IF ( ln_dia_osm .AND. ( iom_use("zdudz_pyc") .OR. iom_use("zdvdz_pyc") ) )                            ln_dia_pyc_shr = .TRUE.
3121      IF ( ln_dia_osm .AND. ( iom_use("zdtdz_pyc") .OR. iom_use("zdsdz_pyc") .OR. iom_use("zdbdz_pyc" ) ) ) ln_dia_pyc_scl = .TRUE.
3122      !
3123      ! Allocate zdfosm arrays
3124      IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
3125      !
3126      IF( ln_osm_mle ) THEN   ! Initialise Fox-Kemper parametrization
3127         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
3128903      IF( ios /= 0 ) CALL ctl_nam( ios, 'namosm_mle in reference namelist' )
3129         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
3130904      IF( ios >  0 ) CALL ctl_nam( ios, 'namosm_mle in configuration namelist' )
3131         IF(lwm) WRITE ( numond, namosm_mle )
3132         !
3133         IF(lwp) THEN   ! Namelist print
3134            WRITE(numout,*)
3135            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
3136            WRITE(numout,*) '~~~~~~~~~~~~~'
3137            WRITE(numout,*) '   Namelist namosm_mle : '
3138            WRITE(numout,*) '      MLE type: =0 standard Fox-Kemper ; =1 new formulation   nn_osm_mle        = ', nn_osm_mle
3139            WRITE(numout,*) '      Magnitude of the MLE (typical value: 0.06 to 0.08)      rn_osm_mle_ce     = ', rn_osm_mle_ce
3140            WRITE(numout,*) '      Scale of ML front (ML radius of deform.) (nn_osm_mle=0) rn_osm_mle_lf     = ', rn_osm_mle_lf,    &
3141               &            'm'
3142            WRITE(numout,*) '      Maximum time scale of MLE                (nn_osm_mle=0) rn_osm_mle_time   = ',   &
3143               &            rn_osm_mle_time, 's'
3144            WRITE(numout,*) '      Reference latitude (deg) of MLE coef.    (nn_osm_mle=1) rn_osm_mle_lat    = ', rn_osm_mle_lat,   &
3145               &            'deg'
3146            WRITE(numout,*) '      Density difference used to define ML for FK             rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c
3147            WRITE(numout,*) '      Threshold used to define MLE for FK                     rn_osm_mle_thresh = ',   &
3148               &            rn_osm_mle_thresh, 'm^2/s'
3149            WRITE(numout,*) '      Timescale for OSM-FK                                    rn_osm_mle_tau    = ', rn_osm_mle_tau, 's'
3150            WRITE(numout,*) '      Switch to limit hmle                                    ln_osm_hmle_limit = ', ln_osm_hmle_limit
3151            WRITE(numout,*) '      hmle limit (fraction of zmld) (ln_osm_hmle_limit = .T.) rn_osm_hmle_limit = ', rn_osm_hmle_limit
3152         END IF
3153      END IF
3154      !
3155      IF(lwp) THEN
3156         WRITE(numout,*)
3157         IF ( ln_osm_mle ) THEN
3158            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
3159            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
3160            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation'
3161         ELSE
3162            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
3163         END IF
3164      END IF
3165      !
3166      IF( ln_osm_mle ) THEN   ! MLE initialisation
3167         !
3168         rb_c = grav * rn_osm_mle_rho_c / rho0   ! Mixed Layer buoyancy criteria
3169         IF(lwp) WRITE(numout,*)
3170         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
3171         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
3172         !
3173         IF( nn_osm_mle == 1 ) THEN
3174            rc_f = rn_osm_mle_ce / ( 5e3_wp * 2.0_wp * omega * SIN( rad * rn_osm_mle_lat ) )
3175         END IF
3176         ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
3177         z1_t2 = 2e-5_wp
3178         DO_2D( 0, 0, 0, 0 )
3179            r1_ft(ji,jj) = MIN( 1.0_wp / ( ABS( ff_t(ji,jj)) + epsln ), ABS( ff_t(ji,jj) ) / z1_t2**2 )
3180         END_2D
3181         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
3182         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
3183         !
3184      END IF
3185      !
3186      CALL osm_rst( nit000, Kmm,  'READ' )   ! Read or initialize hbl, dh, hmle
3187      !
3188      IF ( ln_zdfddm ) THEN
3189         IF(lwp) THEN
3190            WRITE(numout,*)
3191            WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
3192            WRITE(numout,*) '    CAUTION : done in routine zdfosm, not in routine zdfddm '
3193         END IF
3194      END IF
3195      !
3196      ! Set constants not in namelist
3197      ! -----------------------------
3198      IF(lwp) THEN
3199         WRITE(numout,*)
3200      END IF
3201      !
3202      dstokes(:,:) = pp_large
3203      IF (nn_osm_wave == 0) THEN
3204         dstokes(:,:) = rn_osm_dstokes
3205      END IF
3206      !
3207      ! Horizontal average : initialization of weighting arrays
3208      ! -------------------
3209      SELECT CASE ( nn_ave )
3210      CASE ( 0 )                ! no horizontal average
3211         IF(lwp) WRITE(numout,*) '          no horizontal average on avt'
3212         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
3213         ! Weighting mean arrays etmean
3214         !           ( 1  1 )
3215         ! avt = 1/4 ( 1  1 )
3216         !
3217         etmean(:,:,:) = 0.0_wp
3218         !
3219         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
3220            etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.0_wp, umask(ji-1,jj,  jk) + umask(ji,jj,jk) +   &
3221               &                                              vmask(ji,  jj-1,jk) + vmask(ji,jj,jk) )
3222         END_3D
3223      CASE ( 1 )                ! horizontal average
3224         IF(lwp) WRITE(numout,*) '          horizontal average on avt'
3225         ! Weighting mean arrays etmean
3226         !           ( 1/2  1  1/2 )
3227         ! avt = 1/8 ( 1    2  1   )
3228         !           ( 1/2  1  1/2 )
3229         etmean(:,:,:) = 0.0_wp
3230         !
3231         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
3232            etmean(ji,jj,jk) = tmask(ji, jj,jk) / MAX( 1.0_wp, 2.0_wp *   tmask(ji,jj,jk) +                               &
3233               &                                               0.5_wp * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) +     &
3234               &                                                          tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) +   &
3235               &                                               1.0_wp * ( tmask(ji-1,jj,  jk) + tmask(ji,  jj+1,jk) +     &
3236               &                                                          tmask(ji,  jj-1,jk) + tmask(ji+1,jj,  jk) ) )
3237         END_3D
3238      CASE DEFAULT
3239         WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
3240         CALL ctl_stop( ctmp1 )
3241      END SELECT
3242      !
3243      ! Initialization of vertical eddy coef. to the background value
3244      ! -------------------------------------------------------------
3245      DO jk = 1, jpk
3246         avt(:,:,jk) = avtb(jk) * tmask(:,:,jk)
3247      END DO
3248      !
3249      ! Zero the surface flux for non local term and osm mixed layer depth
3250      ! ------------------------------------------------------------------
3251      ghamt(:,:,:) = 0.0_wp
3252      ghams(:,:,:) = 0.0_wp
3253      ghamu(:,:,:) = 0.0_wp
3254      ghamv(:,:,:) = 0.0_wp
3255      !
3256      IF ( ln_dia_osm ) osmdia2d(:,:) = 0.0_wp   ! Initialise auxiliary array for diagnostic output
3257      !
3258      IF( ln_timing ) CALL timing_stop('zdf_osm_init')
3259      !
3260   END SUBROUTINE zdf_osm_init
3261
3262   SUBROUTINE osm_rst( kt, Kmm, cdrw )
3263      !!---------------------------------------------------------------------
3264      !!                   ***  ROUTINE osm_rst  ***
3265      !!
3266      !! ** Purpose :   Read or write BL fields in restart file
3267      !!
3268      !! ** Method  :   use of IOM library. If the restart does not contain
3269      !!                required fields, they are recomputed from stratification
3270      !!
3271      !!----------------------------------------------------------------------
3272      INTEGER         , INTENT(in   ) ::   kt     ! Ocean time step index
3273      INTEGER         , INTENT(in   ) ::   Kmm    ! Ocean time level index (middle)
3274      CHARACTER(len=*), INTENT(in   ) ::   cdrw   ! "READ"/"WRITE" flag
3275      !
3276      ! Local variables
3277      INTEGER  ::   id1, id2, id3                 ! iom enquiry index
3278      INTEGER  ::   ji, jj, jk                    ! Dummy loop indices
3279      INTEGER  ::   iiki, ikt                     ! Local integer
3280      REAL(wp) ::   zhbf                          ! Tempory scalars
3281      REAL(wp) ::   zN2_c                         ! Local scalar
3282      REAL(wp) ::   rho_c = 0.01_wp               ! Density criterion for mixed layer depth
3283      INTEGER, DIMENSION(jpi,jpj) ::   imld_rst   ! Level of mixed-layer depth (pycnocline top)
3284      !
3285      IF( ln_timing ) CALL timing_start('osm_rst')
3286      !
3287      !!-----------------------------------------------------------------------------
3288      ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
3289      !!-----------------------------------------------------------------------------
3290      IF( TRIM(cdrw) == 'READ' .AND. ln_rstart) THEN
3291         id1 = iom_varid( numror, 'wn', ldstop = .FALSE. )
3292         IF( id1 > 0 ) THEN   ! 'wn' exists; read
3293            CALL iom_get( numror, jpdom_auto, 'wn', ww )
3294            WRITE(numout,*) ' ===>>>> :  wn read from restart file'
3295         ELSE
3296            ww(:,:,:) = 0.0_wp
3297            WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
3298         END IF
3299         !
3300         id1 = iom_varid( numror, 'hbl', ldstop = .FALSE. )
3301         id2 = iom_varid( numror, 'dh',  ldstop = .FALSE. )
3302         IF( id1 > 0 .AND. id2 > 0 ) THEN   ! 'hbl' exists; read and return
3303            CALL iom_get( numror, jpdom_auto, 'hbl', hbl  )
3304            CALL iom_get( numror, jpdom_auto, 'dh',  dh   )
3305            hml(:,:) = hbl(:,:) - dh(:,:)   ! Initialise ML depth
3306            WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file'
3307            IF( ln_osm_mle ) THEN
3308               id3 = iom_varid( numror, 'hmle', ldstop = .FALSE. )
3309               IF( id3 > 0 ) THEN
3310                  CALL iom_get( numror, jpdom_auto, 'hmle', hmle )
3311                  WRITE(numout,*) ' ===>>>> :  hmle read from restart file'
3312               ELSE
3313                  WRITE(numout,*) ' ===>>>> :  hmle not found, set to hbl'
3314                  hmle(:,:) = hbl(:,:)   ! Initialise MLE depth
3315               END IF
3316            END IF
3317            RETURN
3318         ELSE   ! 'hbl' & 'dh' not in restart file, recalculate
3319            WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
3320         END IF
3321      END IF
3322      !
3323      !!-----------------------------------------------------------------------------
3324      ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
3325      !!-----------------------------------------------------------------------------
3326      IF ( TRIM(cdrw) == 'WRITE' ) THEN
3327         IF(lwp) WRITE(numout,*) '---- osm-rst ----'
3328         CALL iom_rstput( kt, nitrst, numrow, 'wn',  ww  )
3329         CALL iom_rstput( kt, nitrst, numrow, 'hbl', hbl )
3330         CALL iom_rstput( kt, nitrst, numrow, 'dh',  dh  )
3331         IF ( ln_osm_mle ) THEN
3332            CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle )
3333         END IF
3334         RETURN
3335      END IF
3336      !
3337      !!-----------------------------------------------------------------------------
3338      ! Getting hbl, no restart file with hbl, so calculate from surface stratification
3339      !!-----------------------------------------------------------------------------
3340      IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
3341      ! w-level of the mixing and mixed layers
3342      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )
3343      CALL bn2( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm )
3344      imld_rst(:,:) = nlb10            ! Initialization to the number of w ocean point
3345      hbl(:,:) = 0.0_wp                ! Here hbl used as a dummy variable, integrating vertically N^2
3346      zN2_c = grav * rho_c * r1_rho0   ! Convert density criteria into N^2 criteria
3347      !
3348      hbl(:,:)  = 0.0_wp   ! Here hbl used as a dummy variable, integrating vertically N^2
3349      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
3350         ikt = mbkt(ji,jj)
3351         hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0.0_wp ) * e3w(ji,jj,jk,Kmm)
3352         IF ( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
3353      END_3D
3354      !
3355      DO_2D( 1, 1, 1, 1 )
3356         iiki = MAX( 4, imld_rst(ji,jj) )
3357         hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm  )   ! Turbocline depth
3358         dh(ji,jj)  = e3t(ji,jj,iiki-1,Kmm  )   ! Turbocline depth
3359         hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
3360      END_2D
3361      !
3362      WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
3363      !
3364      IF( ln_osm_mle ) THEN
3365         hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
3366         WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
3367      END IF
3368      !
3369      ww(:,:,:) = 0._wp
3370      WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
3371      !
3372      IF( ln_timing ) CALL timing_stop('osm_rst')
3373      !
3374   END SUBROUTINE osm_rst
3375
3376   SUBROUTINE tra_osm( kt, Kmm, pts, Krhs )
3377      !!----------------------------------------------------------------------
3378      !!                  ***  ROUTINE tra_osm  ***
3379      !!
3380      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
3381      !!
3382      !! ** Method  :   ???
3383      !!
3384      !!----------------------------------------------------------------------
3385      INTEGER                                  , INTENT(in   ) ::   kt          ! Time step index
3386      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs   ! Time level indices
3387      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts         ! Active tracers and RHS of tracer equation
3388      !
3389      ! Local variables
3390      INTEGER                                 ::   ji, jj, jk
3391      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
3392      !
3393      IF( ln_timing ) CALL timing_start('tra_osm')
3394      !
3395      IF ( kt == nit000 ) THEN
3396         IF ( ntile == 0 .OR. ntile == 1 ) THEN   ! Do only on the first tile
3397            IF(lwp) WRITE(numout,*)
3398            IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
3399            IF(lwp) WRITE(numout,*) '~~~~~~~   '
3400         END IF
3401      END IF
3402      !
3403      IF ( l_trdtra ) THEN   ! Save ta and sa trends
3404         ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
3405         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
3406         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
3407      END IF
3408      !
3409      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
3410         pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      &
3411            &                 - (  ghamt(ji,jj,jk  )  &
3412            &                    - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm)
3413         pts(ji,jj,jk,jp_sal,Krhs) =  pts(ji,jj,jk,jp_sal,Krhs)                      &
3414            &                 - (  ghams(ji,jj,jk  )  &
3415            &                    - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
3416      END_3D
3417      !
3418      IF ( l_trdtra ) THEN   ! Save the non-local tracer flux trends for diagnostics
3419         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
3420         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
3421         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt )
3422         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds )
3423         DEALLOCATE( ztrdt, ztrds )
3424      END IF
3425      !
3426      IF ( sn_cfctl%l_prtctl ) THEN
3427         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm  - Ta: ', mask1=tmask,   &
3428            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
3429      END IF
3430      !
3431      IF( ln_timing ) CALL timing_stop('tra_osm')
3432      !
3433   END SUBROUTINE tra_osm
3434
3435   SUBROUTINE trc_osm( kt )   ! Dummy routine
3436      !!----------------------------------------------------------------------
3437      !!                  ***  ROUTINE trc_osm  ***
3438      !!
3439      !! ** Purpose :   compute and add to the passive tracer trend the non-local
3440      !!                 passive tracer flux
3441      !!
3442      !!
3443      !! ** Method  :   ???
3444      !!
3445      !!----------------------------------------------------------------------
3446      INTEGER, INTENT(in) :: kt
3447      !
3448      IF( ln_timing ) CALL timing_start('trc_osm')
3449      !
3450      WRITE(*,*) 'trc_osm: Not written yet', kt
3451      !
3452      IF( ln_timing ) CALL timing_stop('trc_osm')
3453      !
3454   END SUBROUTINE trc_osm
3455
3456   SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs )
3457      !!----------------------------------------------------------------------
3458      !!                  ***  ROUTINE dyn_osm  ***
3459      !!
3460      !! ** Purpose :   compute and add to the velocity trend the non-local flux
3461      !! copied/modified from tra_osm
3462      !!
3463      !! ** Method  :   ???
3464      !!
3465      !!----------------------------------------------------------------------
3466      INTEGER                             , INTENT(in   ) ::   kt          ! Ocean time step index
3467      INTEGER                             , INTENT(in   ) ::   Kmm, Krhs   ! Ocean time level indices
3468      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! Ocean velocities and RHS of momentum equation
3469      !
3470      INTEGER :: ji, jj, jk   ! dummy loop indices
3471      !!----------------------------------------------------------------------
3472      !
3473      IF( ln_timing ) CALL timing_start('dyn_osm')
3474      !
3475      IF ( kt == nit000 ) THEN
3476         IF(lwp) WRITE(numout,*)
3477         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
3478         IF(lwp) WRITE(numout,*) '~~~~~~~   '
3479      END IF
3480      !
3481      ! Code saving tracer trends removed, replace with trdmxl_oce
3482      !
3483      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Add non-local u and v fluxes
3484         puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( ghamu(ji,jj,jk  ) -   &
3485            &                                         ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm)
3486         pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( ghamv(ji,jj,jk  ) -   &
3487            &                                         ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm)
3488      END_3D
3489      !
3490      ! Code for saving tracer trends removed
3491      !
3492      IF( ln_timing ) CALL timing_stop('dyn_osm')
3493      !
3494   END SUBROUTINE dyn_osm
3495
3496   !!======================================================================
3497
3498END MODULE zdfosm
Note: See TracBrowser for help on using the repository browser.