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

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

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/ZDF/zdfosm.F90 @ 14219

Last change on this file since 14219 was 14219, checked in by mcastril, 4 years ago

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