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

Last change on this file since 14555 was 14555, 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: adoption of changeset [14441] (ticket #2353)

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