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

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

source: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/ZDF/zdfosm.F90 @ 13403

Last change on this file since 13403 was 13403, checked in by dancopsey, 4 years ago

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