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

Last change on this file since 14759 was 14759, checked in by agn, 7 months ago

initialiose zlarge in zdf_osm_init

  • Property svn:keywords set to Id
File size: 203.5 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(:,:) = zlarge
2057      zshear(2:jpim1,2:jpjm1) = 0._wp
2058#ifdef key_osm_debug
2059      IF(narea==nn_narea_db) THEN
2060         ji=iloc_db; jj=jloc_db
2061         WRITE(narea+100,'(a,g11.3)') &
2062              & 'zdf_osm_osbl_state start: zekman=', zekman(ji,jj)
2063         FLUSH(narea+100)
2064      END IF
2065#endif
2066      j_ddh(2:jpim1,2:jpjm1) = 1     
2067
2068      DO jj = 2, jpjm1
2069         DO ji = 2, jpim1
2070            IF ( lconv(ji,jj) ) THEN
2071               IF ( zdb_bl(ji,jj) > 0._wp ) THEN
2072                  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 &
2073                       & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp )
2074
2075                  IF ( ff_t(ji,jj) >= 0._wp ) THEN
2076                     !  Northern Hemisphere
2077                     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 )
2078                  ELSE
2079                     !  Southern Hemisphere
2080                     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 )
2081                  ENDIF
2082                  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 ) )
2083#ifdef key_osm_debug
2084                  ! IF(narea==nn_narea_db)THEN
2085                  !    WRITE(narea+100,'(2(a,i10.4))')'ji',ji,'jj',jj
2086                  !    WRITE(narea+100,'(2(a,i10.4))')'iloc_db',iloc_db,'jloc_db',jloc_db
2087                  !    WRITE(narea+100,'(2(a,i10.4))')'iloc_db+',mi0(nn_idb),'jloc_db+',mj0(nn_jdb)
2088                  !    FLUSH(narea+100)
2089                  ! END IF
2090                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2091                     WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear=',zshear(ji,jj)
2092                     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)
2093                     FLUSH(narea+100)
2094                  END IF
2095#endif
2096                  ! Stability Dependence
2097                  zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75 * MAX(0._wp,( zri_b(ji,jj) - zri_c ) / zri_c ) )
2098#ifdef key_osm_debug
2099                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2100                     WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear inc ri part=',zshear(ji,jj)
2101                     FLUSH(narea+100)
2102                  END IF
2103#endif
2104
2105!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2106                  ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when  !
2107                  ! full code available                                          !
2108!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2109                  IF ( zshear(ji,jj) > 1.e-10 ) THEN
2110                     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
2111                        ! Growing shear layer
2112                        j_ddh(ji,jj) = 0
2113                        lshear(ji,jj) = .TRUE.
2114                     ELSE
2115                        j_ddh(ji,jj) = 1
2116                        !                IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN
2117                        ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer.
2118                        lshear(ji,jj) = .TRUE.
2119                        !                ELSE
2120                     ENDIF
2121                  ELSE
2122                     j_ddh(ji,jj) = 2
2123                     lshear(ji,jj) = .FALSE.
2124                  ENDIF
2125                  ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline.
2126                  !                  zshear(ji,jj) = 0.5 * zshear(ji,jj)
2127                  !                  lshear(ji,jj) = .FALSE.
2128                  !                ENDIF               
2129               ELSE                ! zdb_bl test, note zshear set to zero
2130                  j_ddh(ji,jj) = 2
2131                  lshear(ji,jj) = .FALSE.
2132               ENDIF
2133            ENDIF
2134         END DO
2135      END DO
2136
2137      ! Calculate entrainment buoyancy flux due to surface fluxes.
2138
2139      DO jj = 2, jpjm1
2140         DO ji = 2, jpim1
2141            IF ( lconv(ji,jj) ) THEN
2142               zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
2143               zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
2144               zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
2145               zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
2146               IF (nn_osm_SD_reduce > 0 ) THEN
2147                  ! Effective Stokes drift already reduced from surface value
2148                  zr_stokes = 1.0_wp
2149               ELSE
2150                  ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
2151                  ! requires further reduction where BL is deep
2152                  zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
2153                       &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
2154               END IF
2155               zwb_ent(ji,jj) = - 2.0 * zalpha_c * zrf_conv * zwbav(ji,jj) &
2156                    &                  - zalpha_s * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
2157                    &         + zr_stokes * ( zalpha_s * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
2158                    &                                         - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
2159               !
2160#ifdef key_osm_debug
2161               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2162                  WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state conv+shear0/lang: zwb_ent=',zwb_ent(ji,jj)
2163                  FLUSH(narea+100)
2164               END IF
2165#endif
2166            ENDIF
2167         END DO ! ji loop
2168      END DO   ! jj loop
2169
2170      zwb_min(:,:) = zlarge
2171
2172      DO jj = 2, jpjm1
2173         DO ji = 2, jpim1
2174            IF ( lshear(ji,jj) ) THEN
2175               IF ( lconv(ji,jj) ) THEN
2176                  ! Unstable OSBL
2177                  zwb_shr = -za_wb_s * zri_b(ji,jj) * zshear(ji,jj)
2178#ifdef key_osm_debug
2179                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2180                     WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zwb_shr: zwb_shr=',zwb_shr
2181                     FLUSH(narea+100)
2182                  END IF
2183#endif
2184                  IF ( j_ddh(ji,jj) == 0 ) THEN
2185
2186                     ! ! Developing shear layer, additional shear production possible.
2187
2188                     !                 zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) /  zhbl(ji,jj), 0._wp )
2189                     !                 zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 )
2190                     !                 zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )
2191
2192                     !                 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 )
2193                     !                 zwb_shr = MAX( zwb_shr, -0.25 * zshear_u )
2194#ifdef key_osm_debug
2195                     IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2196                        WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state j_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, &
2197                             & '  zshear=',zshear(ji,jj),'  zshear_u=', zshear_u
2198                        FLUSH(narea+100)
2199                     END IF
2200#endif
2201
2202                  ENDIF
2203                  zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
2204                  !              zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
2205               ELSE    ! IF ( lconv ) THEN - ENDIF
2206                  ! Stable OSBL  - shear production not coded for first attempt.           
2207               ENDIF  ! lconv
2208            ENDIF  ! lshear
2209            IF ( lconv(ji,jj) ) THEN
2210               ! Unstable OSBL
2211               zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2._wp * zwbav(ji,jj)
2212            ENDIF  ! lconv
2213#ifdef key_osm_debug
2214            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2215               WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',zwb_ent(ji,jj), &
2216                    & '  zwb_min=',zwb_min(ji,jj), '  zwb0tot=', zwb0tot(ji,jj), '  zwbav= ', zwbav(ji,jj)
2217               FLUSH(narea+100)
2218            END IF
2219#endif
2220         END DO   ! ji
2221      END DO     ! jj
2222    END SUBROUTINE zdf_osm_osbl_state
2223
2224
2225    SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv )
2226      !!---------------------------------------------------------------------
2227      !!                   ***  ROUTINE zdf_vertical_average  ***
2228      !!
2229      !! ** Purpose : Determines vertical averages from surface to jnlev.
2230      !!
2231      !! ** Method  : Averages are calculated from the surface to jnlev.
2232      !!              The external level used to calculate differences is ibld+ibld_ext
2233      !!
2234      !!----------------------------------------------------------------------
2235
2236      INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over.
2237      INTEGER, DIMENSION(jpi,jpj) :: jp_ext
2238
2239      ! Alan: do we need zb?
2240      REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb        ! Average temperature and salinity
2241      REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components
2242      REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
2243      REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components.
2244
2245      INTEGER :: jk, ji, jj, ibld_ext
2246      REAL(wp) :: zthick, zthermal, zbeta
2247
2248
2249      zt(2:jpim1,2:jpjm1)  = 0._wp
2250      zs(2:jpim1,2:jpjm1)  = 0._wp
2251      zu(2:jpim1,2:jpjm1)  = 0._wp
2252      zv(2:jpim1,2:jpjm1)  = 0._wp
2253      DO jj = 2, jpjm1                                 !  Vertical slab
2254         DO ji = 2, jpim1
2255            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
2256            zbeta    = rab_n(ji,jj,1,jp_sal)
2257            ! average over depth of boundary layer
2258            zthick = epsln
2259            DO jk = 2, jnlev_av(ji,jj)
2260               zthick = zthick + e3t_n(ji,jj,jk)
2261               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
2262               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
2263               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) &
2264                    &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
2265                    &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
2266               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) &
2267                    &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
2268                    &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
2269            END DO
2270            zt(ji,jj) = zt(ji,jj) / zthick
2271            zs(ji,jj) = zs(ji,jj) / zthick
2272            zu(ji,jj) = zu(ji,jj) / zthick
2273            zv(ji,jj) = zv(ji,jj) / zthick
2274            zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj)
2275            ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj)
2276            IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN      ! ag 09/03
2277 ! Two external levels are available.
2278               zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem)
2279               zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal)
2280               zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) &
2281                    &    / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) )
2282               zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) &
2283                    &   / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) )
2284               zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
2285            ELSE
2286               zdt(ji,jj) = 0._wp
2287               zds(ji,jj) = 0._wp
2288               zdu(ji,jj) = 0._wp
2289               zdv(ji,jj) = 0._wp
2290               zdb(ji,jj) = 0._wp
2291            ENDIF
2292         END DO
2293      END DO
2294    END SUBROUTINE zdf_osm_vertical_average
2295
2296    SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
2297      !!---------------------------------------------------------------------
2298      !!                   ***  ROUTINE zdf_velocity_rotation  ***
2299      !!
2300      !! ** Purpose : Rotates frame of reference of averaged velocity components.
2301      !!
2302      !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w
2303      !!
2304      !!----------------------------------------------------------------------
2305
2306      REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle
2307      REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current
2308      REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline
2309
2310      INTEGER :: ji, jj
2311      REAL(wp) :: ztemp
2312
2313      DO jj = 2, jpjm1
2314         DO ji = 2, jpim1
2315            ztemp = zu(ji,jj)
2316            zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
2317            zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
2318            ztemp = zdu(ji,jj)
2319            zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
2320            zdv(ji,jj) = zdv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
2321         END DO
2322      END DO
2323    END SUBROUTINE zdf_osm_velocity_rotation
2324
2325    SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
2326      !!---------------------------------------------------------------------
2327      !!                   ***  ROUTINE zdf_osm_osbl_state_fk  ***
2328      !!
2329      !! ** 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.
2330      !!  lpyc :: determines whether pycnocline flux-grad relationship needs to be determined
2331      !!  lflux :: determines whether effects of surface flux extend below the base of the OSBL
2332      !!  lmle  :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl.
2333      !!
2334      !! ** Method  :
2335      !!
2336      !!
2337      !!----------------------------------------------------------------------
2338
2339      ! Outputs
2340      LOGICAL,  DIMENSION(jpi,jpj)  :: lpyc, lflux, lmle
2341      REAL(wp), DIMENSION(jpi,jpj)  :: zwb_fk
2342      !
2343      REAL(wp), DIMENSION(jpi,jpj)  :: znd_param
2344      REAL(wp)                      :: zbuoy, ztmp, zpe_mle_layer
2345      REAL(wp)                      :: zpe_mle_ref, zdbdz_mle_int
2346
2347      znd_param(2:jpim1,2:jpjm1) = 0._wp
2348
2349      DO jj = 2, jpjm1
2350         DO ji = 2, jpim1         
2351            ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
2352            zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj)
2353         END DO
2354      END DO
2355      DO jj = 2, jpjm1
2356         DO ji = 2, jpim1
2357            !
2358            IF ( lconv(ji,jj) ) THEN
2359               IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
2360                  zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2361                  zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2362                  zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2363                  zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2364                  ! Calculate potential energies of actual profile and reference profile.
2365                  zpe_mle_layer = 0._wp
2366                  zpe_mle_ref = 0._wp
2367                  zthermal = rab_n(ji,jj,1,jp_tem)
2368                  zbeta    = rab_n(ji,jj,1,jp_sal)
2369
2370                  DO jk = ibld(ji,jj), mld_prof(ji,jj)
2371                     zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) )
2372                     zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
2373                     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)
2374                  END DO
2375                  ! Non-dimensional parameter to diagnose the presence of thermocline
2376
2377                  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) )
2378               ENDIF
2379            ENDIF
2380#ifdef key_osm_debug
2381            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2382               WRITE(narea+100,'(4(a,g11.3))')'start of zdf_osm_osbl_state_fk: zwb_fk=',zwb_fk(ji,jj), &
2383                    & '  znd_param=',znd_param(ji,jj), ' zpe_mle_ref=', zpe_mle_ref,  ' zpe_mle_layer=', zpe_mle_layer
2384               FLUSH(narea+100)
2385            END IF
2386#endif
2387         END DO
2388      END DO
2389
2390      ! Diagnosis
2391      DO jj = 2, jpjm1
2392         DO ji = 2, jpim1
2393            IF ( lconv(ji,jj) ) THEN
2394               IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN
2395                  IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
2396                     ! MLE layer growing
2397                     IF ( znd_param (ji,jj) > 100. ) THEN
2398                        ! Thermocline present
2399                        lflux(ji,jj) = .FALSE.
2400                        lmle(ji,jj) =.FALSE.
2401                     ELSE
2402                        ! Thermocline not present
2403                        lflux(ji,jj) = .TRUE.
2404                        lmle(ji,jj) = .TRUE.
2405                     ENDIF  ! znd_param > 100
2406                     !
2407                     IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
2408                        lpyc(ji,jj) = .FALSE.
2409                     ELSE
2410                        lpyc(ji,jj) = .TRUE.
2411                     ENDIF
2412                  ELSE
2413                     ! MLE layer restricted to OSBL or just below.
2414                     IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
2415                        ! Weak stratification MLE layer can grow.
2416                        lpyc(ji,jj) = .FALSE.
2417                        lflux(ji,jj) = .TRUE.
2418                        lmle(ji,jj) = .TRUE.
2419                     ELSE
2420                        ! Strong stratification
2421                        lpyc(ji,jj) = .TRUE.
2422                        lflux(ji,jj) = .FALSE.
2423                        lmle(ji,jj) = .FALSE.
2424                     ENDIF ! zdb_bl < rn_mle_thresh_bl and
2425                  ENDIF  ! zhmle > 1.2 zhbl
2426               ELSE
2427                  lpyc(ji,jj) = .TRUE.
2428                  lflux(ji,jj) = .FALSE.
2429                  lmle(ji,jj) = .FALSE.
2430                  IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
2431               ENDIF !  -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5
2432            ELSE
2433               ! Stable Boundary Layer
2434               lpyc(ji,jj) = .FALSE.
2435               lflux(ji,jj) = .FALSE.
2436               lmle(ji,jj) = .FALSE.
2437            ENDIF  ! lconv
2438#ifdef key_osm_debug
2439            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2440               WRITE(narea+100,'(3(a,g11.3),/,4(a,l2))')'end of zdf_osm_osbl_state_fk:zwb_ent=',zwb_ent(ji,jj), &
2441                    & '  zhmle=',zhmle(ji,jj), ' zhbl=', zhbl(ji,jj), &
2442                    & ' lpyc= ', lpyc(ji,jj), ' lflux= ', lflux(ji,jj),  ' lmle= ', lmle(ji,jj), ' lconv= ', lconv(ji,jj)
2443               FLUSH(narea+100)
2444            END IF
2445#endif
2446         END DO
2447      END DO
2448    END SUBROUTINE zdf_osm_osbl_state_fk
2449
2450    SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz )
2451      !!---------------------------------------------------------------------
2452      !!                   ***  ROUTINE zdf_osm_external_gradients  ***
2453      !!
2454      !! ** Purpose : Calculates the gradients below the OSBL
2455      !!
2456      !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient.
2457      !!
2458      !!----------------------------------------------------------------------
2459
2460      INTEGER, DIMENSION(jpi,jpj)  :: jbase
2461      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy.
2462
2463      INTEGER :: jj, ji, jkb, jkb1
2464      REAL(wp) :: zthermal, zbeta
2465
2466
2467      DO jj = 2, jpjm1
2468         DO ji = 2, jpim1
2469            IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
2470               zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
2471               zbeta    = rab_n(ji,jj,1,jp_sal)
2472               jkb = jbase(ji,jj)
2473               jkb1 = MIN(jkb + 1, mbkt(ji,jj))
2474               zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
2475                    &   / e3w_n(ji,jj,jkb1)
2476               zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
2477                    &   / e3w_n(ji,jj,jkb1)
2478               zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
2479            ELSE
2480               zdtdz(ji,jj) = 0._wp
2481               zdsdz(ji,jj) = 0._wp
2482               zdbdz(ji,jj) = 0._wp
2483            END IF
2484         END DO
2485      END DO
2486    END SUBROUTINE zdf_osm_external_gradients
2487
2488    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha )
2489
2490      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline
2491      REAL(wp), DIMENSION(jpi,jpj) :: zalpha
2492
2493      INTEGER :: jk, jj, ji
2494      REAL(wp) :: ztgrad, zsgrad, zbgrad
2495      REAL(wp) :: zgamma_b_nd, znd
2496      REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc
2497      REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
2498
2499      DO jj = 2, jpjm1
2500         DO ji = 2, jpim1
2501            IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
2502               IF ( lconv(ji,jj) ) THEN  ! convective conditions
2503                  IF ( lpyc(ji,jj) ) THEN
2504                     zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
2505                     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 ) )
2506                     zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp )
2507
2508                     ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
2509!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2510                     ! Commented lines in this section are not needed in new code, once tested !
2511                     ! can be removed                                                          !
2512!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2513                     !                   ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
2514                     !                   zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
2515                     zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)
2516                     zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
2517                     DO jk = 2, ibld(ji,jj)
2518                        znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp
2519                        IF ( znd <= zzeta_m ) THEN
2520                           !                        zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &
2521                           !                &        EXP( -6.0 * ( znd -zzeta_m )**2 )
2522                           !                        zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &
2523                           !                                  & EXP( -6.0 * ( znd -zzeta_m )**2 )
2524                           zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &
2525                                & EXP( -6.0 * ( znd -zzeta_m )**2 )
2526                        ELSE
2527                           !                         zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
2528                           !                         zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
2529                           zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
2530                        ENDIF
2531                     END DO
2532#ifdef key_osm_debug
2533                     IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2534                        WRITE(narea+100,'(a,/,3(a,g11.3),/,2(a,g11.3),/)')'end of zdf_osm_pycnocline_scalar_profiles:lconv=lpyc=T',&
2535                             & 'zzeta_m=', zzeta_m, ' zalpha=', zalpha(ji,jj), ' ztmp=', ztmp,&
2536                             & ' zbgrad=', zbgrad, ' zgamma_b_nd=', zgamma_b_nd
2537                        FLUSH(narea+100)
2538                     END IF
2539#endif
2540                  ENDIF ! if no pycnocline pycnocline gradients set to zero
2541               ELSE
2542                  ! stable conditions
2543                  ! if pycnocline profile only defined when depth steady of increasing.
2544                  IF ( zdhdt(ji,jj) > 0.0 ) THEN        ! Depth increasing, or steady.
2545                     IF ( zdb_bl(ji,jj) > 0._wp ) THEN
2546                        IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
2547                           ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
2548                           ztgrad = zdt_bl(ji,jj) * ztmp
2549                           zsgrad = zds_bl(ji,jj) * ztmp
2550                           zbgrad = zdb_bl(ji,jj) * ztmp
2551                           DO jk = 2, ibld(ji,jj)
2552                              znd = gdepw_n(ji,jj,jk) * ztmp
2553                              zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
2554                              zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
2555                              zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
2556                           END DO
2557                        ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
2558                           ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
2559                           ztgrad = zdt_bl(ji,jj) * ztmp
2560                           zsgrad = zds_bl(ji,jj) * ztmp
2561                           zbgrad = zdb_bl(ji,jj) * ztmp
2562                           DO jk = 2, ibld(ji,jj)
2563                              znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp
2564                              zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
2565                              zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
2566                              zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
2567                           END DO
2568                        ENDIF ! IF (zhol >=0.5)
2569#ifdef key_osm_debug
2570                        IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2571                           WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_pycnocline_scalar_profiles:lconv=F ztgrad=',&
2572                                & ztgrad, ' zsgrad=', zsgrad, ' zbgrad=', zbgrad
2573                           FLUSH(narea+100)
2574                        END IF
2575#endif
2576                     ENDIF    ! IF (zdb_bl> 0.)
2577                  ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
2578               ENDIF          ! IF (lconv)
2579            ENDIF      ! IF ( ibld(ji,jj) < mbkt(ji,jj) )
2580         END DO
2581      END DO
2582
2583    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
2584
2585    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
2586      !!---------------------------------------------------------------------
2587      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  ***
2588      !!
2589      !! ** Purpose : Calculates velocity shear in the pycnocline
2590      !!
2591      !! ** Method  :
2592      !!
2593      !!----------------------------------------------------------------------
2594
2595      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
2596
2597      INTEGER :: jk, jj, ji
2598      REAL(wp) :: zugrad, zvgrad, znd
2599      REAL(wp) :: zzeta_v = 0.45
2600      !
2601      DO jj = 2, jpjm1
2602         DO ji = 2, jpim1
2603            !
2604            IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
2605               IF ( lconv (ji,jj) ) THEN
2606                  ! Unstable conditions. Shouldn;t be needed with no pycnocline code.
2607                  !                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
2608                  !                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
2609                  !                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
2610                  !Alan is this right?
2611                  !                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
2612                  !                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
2613                  !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
2614                  !                       &      )/ (zdh(ji,jj)  + epsln )
2615                  !                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
2616                  !                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
2617                  !                     IF ( znd <= 0.0 ) THEN
2618                  !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
2619                  !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
2620                  !                     ELSE
2621                  !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
2622                  !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
2623                  !                     ENDIF
2624                  !                  END DO
2625               ELSE
2626                  ! stable conditions
2627                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
2628                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
2629                  DO jk = 2, ibld(ji,jj)
2630                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
2631                     IF ( znd < 1.0 ) THEN
2632                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
2633                     ELSE
2634                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
2635                     ENDIF
2636                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
2637                  END DO
2638               ENDIF
2639               !
2640            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
2641         END DO
2642      END DO
2643    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
2644
2645    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt )
2646      !!---------------------------------------------------------------------
2647      !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
2648      !!
2649      !! ** Purpose : Calculates the rate at which hbl changes.
2650      !!
2651      !! ** Method  :
2652      !!
2653      !!----------------------------------------------------------------------
2654
2655      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt        ! Rate of change of hbl
2656
2657      INTEGER :: jj, ji
2658      REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi
2659      REAL(wp) :: zvel_max,zddhdt
2660      REAL(wp), PARAMETER :: zzeta_m = 0.3
2661      REAL(wp), PARAMETER :: zgamma_c = 2.0
2662      REAL(wp), PARAMETER :: zdhoh = 0.1
2663      REAL(wp), PARAMETER :: zalpha_b = 0.3
2664      REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
2665
2666      DO jj = 2, jpjm1
2667         DO ji = 2, jpim1
2668
2669            IF ( lshear(ji,jj) ) THEN
2670               IF ( lconv(ji,jj) ) THEN    ! Convective
2671
2672                  IF ( ln_osm_mle ) THEN
2673
2674                     IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
2675                        ! Fox-Kemper buoyancy flux average over OSBL
2676                        zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  &
2677                             (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
2678                     ELSE
2679                        zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
2680                     ENDIF
2681                     zvel_max = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
2682                     IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
2683                        ! OSBL is deepening, entrainment > restratification
2684                        IF ( zdb_bl(ji,jj) > 1.0e-15 ) THEN
2685                           zgamma_b_nd = MAX( zdbdz_bl_ext(ji,jj), 0._wp ) * zdh(ji,jj) /  ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
2686                           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)
2687                           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 )
2688                           zpsi = zalpha_b * MAX ( zpsi, 0._wp )
2689                           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 ) )
2690#ifdef key_osm_debug
2691                           IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2692                              WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt, OSBL is deepening, entrainment > restratification:  zdhdt=',zdhdt(ji,jj)
2693                              WRITE(narea+100,'(3(a,g11.3))') '  zpsi=',zpsi, '  zgamma_b_nd=', zgamma_b_nd, '  zdh=', zdh(ji,jj)
2694                              FLUSH(narea+100)
2695                           END IF
2696#endif
2697                           IF ( j_ddh(ji,jj) == 1 ) THEN
2698                              IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
2699                                 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.2