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/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF – NEMO

source: NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90 @ 14519

Last change on this file since 14519 was 14519, checked in by agn, 3 years ago

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