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

Last change on this file since 14564 was 14564, checked in by smueller, 7 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 [14410,14411] (ticket #2353)

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