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

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

Extra diagnostics to elucidate pyc issues

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