source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90 @ 13859

Last change on this file since 13859 was 13859, checked in by agn, 6 months ago

Added Alan's commit c46eae7f7fd788 + corr to zznd+pyc def as email 1613 23 Nov 20

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