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/trunk/src/OCE/ZDF – NEMO

source: NEMO/trunk/src/OCE/ZDF/zdfosm.F90

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

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