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

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

code uses mi0/mi1 & mj0/mj1 tospecify nn_area

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