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_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90 @ 14149

Last change on this file since 14149 was 14149, checked in by smueller, 3 years ago

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