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

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

first try of Alan's new code

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