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/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF – NEMO

source: NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90 @ 14729

Last change on this file since 14729 was 14729, checked in by smueller, 3 years ago

Synchronisation of the OSMOSIS boundary layer scheme with the version developed in branch /NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0: transfer of changesets [14645,14646] (ticket #2353)

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