source: NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfosm.F90 @ 12266

Last change on this file since 12266 was 12266, checked in by cguiavarch, 11 months ago

Debug

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