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

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

Simplify undulation term

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