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

Last change on this file since 14571 was 14571, checked in by smueller, 5 months ago

Synchronisation of the OSMOSIS boundary layer scheme with the version developed in branch /NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0: transfer of changesets [14406,14518,14521,14534,14539,14540] (ticket #2353)

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