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/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF – NEMO

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90 @ 14636

Last change on this file since 14636 was 14636, checked in by hadcv, 3 years ago

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