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 @ 14280

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

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