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

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

Ensure pb_coup is defined and used when lcoup=T; include factor of 0.5 in zmsk[uv]

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