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

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

source: NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90 @ 14305

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

Reduction of memory usage by subroutine zdf_osm of module zdfosm (ticket #2353)

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