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

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

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

Last change on this file since 14405 was 14405, checked in by agn, 3 years ago

Alan's new code for when shear is switched on

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