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

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

seems to be minimal zeroing required

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