New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfosm.F90 in NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF – NEMO

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

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

Correct loop in tramsport term in flux-gradient relationship so starts at ji,jj = 2,2

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