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/UKMO/NEMO_4.0.2_FKOSM/src/OCE/ZDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/ZDF/zdfosm.F90

Last change on this file was 12912, checked in by cguiavarch, 4 years ago

OSMOSIS and Fox-Kemper changes (merged from NEMO_4.0.1_FKOSM)

File size: 127.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_init  : initialization, namelist read, and parameters control
43   !!   osm_rst       : read (or initialize) and write osmosis restart fields
44   !!   tra_osm       : compute and add to the T & S trend the non-local flux
45   !!   trc_osm       : compute and add to the passive tracer trend the non-local flux (TBD)
46   !!   dyn_osm       : compute and add to u & v trensd the non-local flux
47   !!
48   !! Subroutines in revised code.
49   !!----------------------------------------------------------------------
50   USE oce            ! ocean dynamics and active tracers
51                      ! uses wn from previous time step (which is now wb) to calculate hbl
52   USE dom_oce        ! ocean space and time domain
53   USE zdf_oce        ! ocean vertical physics
54   USE sbc_oce        ! surface boundary condition: ocean
55   USE sbcwave        ! surface wave parameters
56   USE phycst         ! physical constants
57   USE eosbn2         ! equation of state
58   USE traqsr         ! details of solar radiation absorption
59   USE zdfddm         ! double diffusion mixing (avs array)
60   USE iom            ! I/O library
61   USE lib_mpp        ! MPP library
62   USE trd_oce        ! ocean trends definition
63   USE trdtra         ! tracers trends
64   !
65   USE in_out_manager ! I/O manager
66   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
67   USE prtctl         ! Print control
68   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
69
70   IMPLICIT NONE
71   PRIVATE
72
73   PUBLIC   zdf_osm       ! routine called by step.F90
74   PUBLIC   zdf_osm_init  ! routine called by nemogcm.F90
75   PUBLIC   osm_rst       ! routine called by step.F90
76   PUBLIC   tra_osm       ! routine called by step.F90
77   PUBLIC   trc_osm       ! routine called by trcstp.F90
78   PUBLIC   dyn_osm       ! routine called by step.F90
79
80   PUBLIC   ln_osm_mle    ! logical needed by tra_mle_init in tramle.F90
81
82   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux
83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv    !: non-local v-momentum flux
84   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamt    !: non-local temperature flux (gamma/<ws>o)
85   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghams    !: non-local salinity flux (gamma/<ws>o)
86   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt
87   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth
88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! depth of pycnocline
89   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml      ! ML depth
90   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift.
91
92   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   r1_ft    ! inverse of the modified Coriolis parameter at t-pts
93   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle     ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization
94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle ! zonal buoyancy gradient in ML
95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle ! meridional buoyancy gradient in ML
96   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof ! level of base of MLE layer.
97
98   !                      !!** Namelist  namzdf_osm  **
99   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la
100
101   LOGICAL  ::   ln_osm_mle           !: flag to activate the Mixed Layer Eddy (MLE) parameterisation
102
103   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number
104   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift
105   REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs
106   INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt
107   INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
108   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la
109
110
111   LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing
112   REAL(wp) ::   rn_riinfty   = 0.7     ! local Richardson Number limit for shear instability
113   REAL(wp) ::   rn_difri    =  0.005   ! maximum shear mixing at Rig = 0    (m2/s)
114   LOGICAL  ::   ln_convmix  = .true.   ! Convective instability mixing
115   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s)
116
117! OSMOSIS mixed layer eddy parametrization constants
118   INTEGER  ::   nn_osm_mle             ! = 0/1 flag for horizontal average on avt
119   REAL(wp) ::   rn_osm_mle_ce           ! MLE coefficient
120   !                                        ! parameters used in nn_osm_mle = 0 case
121   REAL(wp) ::   rn_osm_mle_lf               ! typical scale of mixed layer front
122   REAL(wp) ::   rn_osm_mle_time             ! time scale for mixing momentum across the mixed layer
123   !                                        ! parameters used in nn_osm_mle = 1 case
124   REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front
125   REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK
126   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation
127   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld
128   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case
129   REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base.
130   REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE.
131
132
133   !                                    !!! ** General constants  **
134   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero
135   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw
136   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3
137   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3
138
139   INTEGER :: idebug = 236
140   INTEGER :: jdebug = 228
141   !!----------------------------------------------------------------------
142   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
143   !! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $
144   !! Software governed by the CeCILL license (see ./LICENSE)
145   !!----------------------------------------------------------------------
146CONTAINS
147
148   INTEGER FUNCTION zdf_osm_alloc()
149      !!----------------------------------------------------------------------
150      !!                 ***  FUNCTION zdf_osm_alloc  ***
151      !!----------------------------------------------------------------------
152     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &
153          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
154          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
155
156     ALLOCATE(  hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), &
157          &       mld_prof(jpi,jpj), STAT= zdf_osm_alloc )
158
159!     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &    ! would ths be better ?
160!          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
161!          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
162!     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
163!     IF ( ln_osm_mle ) THEN
164!        Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc )
165!        IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays')
166!     ENDIF
167
168     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
169     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
170   END FUNCTION zdf_osm_alloc
171
172
173   SUBROUTINE zdf_osm( kt, p_avm, p_avt )
174      !!----------------------------------------------------------------------
175      !!                   ***  ROUTINE zdf_osm  ***
176      !!
177      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
178      !!      coefficients and non local mixing using the OSMOSIS scheme
179      !!
180      !! ** Method :   The boundary layer depth hosm is diagnosed at tracer points
181      !!      from profiles of buoyancy, and shear, and the surface forcing.
182      !!      Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
183      !!
184      !!                      Kx =  hosm  Wx(sigma) G(sigma)
185      !!
186      !!             and the non local term ghamt = Cs / Ws(sigma) / hosm
187      !!      Below hosm  the coefficients are the sum of mixing due to internal waves
188      !!      shear instability and double diffusion.
189      !!
190      !!      -1- Compute the now interior vertical mixing coefficients at all depths.
191      !!      -2- Diagnose the boundary layer depth.
192      !!      -3- Compute the now boundary layer vertical mixing coefficients.
193      !!      -4- Compute the now vertical eddy vicosity and diffusivity.
194      !!      -5- Smoothing
195      !!
196      !!        N.B. The computation is done from jk=2 to jpkm1
197      !!             Surface value of avt are set once a time to zero
198      !!             in routine zdf_osm_init.
199      !!
200      !! ** Action  :   update the non-local terms ghamts
201      !!                update avt (before vertical eddy coef.)
202      !!
203      !! References : Large W.G., Mc Williams J.C. and Doney S.C.
204      !!         Reviews of Geophysics, 32, 4, November 1994
205      !!         Comments in the code refer to this paper, particularly
206      !!         the equation number. (LMD94, here after)
207      !!----------------------------------------------------------------------
208      INTEGER                   , INTENT(in   ) ::   kt            ! ocean time step
209      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points)
210      !!
211      INTEGER ::   ji, jj, jk                   ! dummy loop indices
212
213      INTEGER ::   jl                   ! dummy loop indices
214
215      INTEGER ::   ikbot, jkmax, jkm1, jkp2     !
216
217      REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube      !
218      REAL(wp) ::   zbeta, zthermal                                  !
219      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
220      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      !
221
222      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density
223      INTEGER  ::   jm                          ! dummy loop indices
224      REAL(wp) ::   zr1, zr2, zr3, zr4, zrhop   ! Compression terms
225      REAL(wp) ::   zflag, zrn2, zdep21, zdep32, zdep43
226      REAL(wp) ::   zesh2, zri, zfri            ! Interior richardson mixing
227      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
228      REAL(wp) :: zt,zs,zu,zv,zrh               ! variables used in constructing averages
229! Scales
230      REAL(wp), DIMENSION(jpi,jpj) :: zrad0     ! Surface solar temperature flux (deg m/s)
231      REAL(wp), DIMENSION(jpi,jpj) :: zradh     ! Radiative flux at bl base (Buoyancy units)
232      REAL(wp), DIMENSION(jpi,jpj) :: zradav    ! Radiative flux, bl average (Buoyancy Units)
233      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
234      REAL(wp), DIMENSION(jpi,jpj) :: zwstrl    ! Langmuir velocity scale
235      REAL(wp), DIMENSION(jpi,jpj) :: zvstr     ! Velocity scale that ends to zustar for large Langmuir number.
236      REAL(wp), DIMENSION(jpi,jpj) :: zwstrc    ! Convective velocity scale
237      REAL(wp), DIMENSION(jpi,jpj) :: zuw0      ! Surface u-momentum flux
238      REAL(wp), DIMENSION(jpi,jpj) :: zvw0      ! Surface v-momentum flux
239      REAL(wp), DIMENSION(jpi,jpj) :: zwth0     ! Surface heat flux (Kinematic)
240      REAL(wp), DIMENSION(jpi,jpj) :: zws0      ! Surface freshwater flux
241      REAL(wp), DIMENSION(jpi,jpj) :: zwb0      ! Surface buoyancy flux
242      REAL(wp), DIMENSION(jpi,jpj) :: zwthav    ! Heat flux - bl average
243      REAL(wp), DIMENSION(jpi,jpj) :: zwsav     ! freshwater flux - bl average
244      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average
245      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux
246
247
248      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b  ! MLE buoyancy flux averaged over OSBL
249      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk    ! max MLE buoyancy flux
250      REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff
251      REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle  ! velocity scale for dhdt with stable ML and FK
252
253      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift
254      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number
255      REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
256      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
257      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer
258      LOGICAL, DIMENSION(jpi,jpj)  :: lconv     ! unstable/stable bl
259
260      ! mixed-layer variables
261
262      INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
263      INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
264
265      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
266      REAL(wp) :: zugrad,zvgrad        ! temporary variables for calculating pycnocline shear
267
268      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid
269      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid
270
271      REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid
272      REAL(wp), DIMENSION(jpi,jpj) :: zmld  ! ML depth on grid
273
274      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid
275      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
276      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2                                    ! correction to dhdt due to internal structure.
277      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext              ! external temperature/salinity and buoyancy gradients
278
279      REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy      ! horizontal gradients for Fox-Kemper parametrization.
280
281      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer
282      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer
283      REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer
284      REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
285      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
286      REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline
287      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline
288      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline
289      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline
290      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc    ! u-shear across the pycnocline
291      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc    ! v-shear across the pycnocline
292
293      ! Flux-gradient relationship variables
294
295      REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale.
296
297      REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity
298      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
299      REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms.
300      REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
301
302      ! For calculating Ri#-dependent mixing
303      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du   ! u-shear^2
304      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv   ! v-shear^2
305      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
306
307      ! Temporary variables
308      INTEGER :: inhml
309      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
310      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables
311      REAL(wp) :: zthick, zz0, zz1 ! temporary variables
312      REAL(wp) :: zvel_max, zhbl_s ! temporary variables
313      REAL(wp) :: zfac             ! temporary variable
314      REAL(wp) :: zus_x, zus_y     ! temporary Stokes drift
315      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
316      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
317
318      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme
319      REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt
320      REAL(wp) :: zzeta_s = 0._wp
321      REAL(wp) :: zzeta_v = 0.46
322      REAL(wp) :: zabsstke
323
324      ! For debugging
325      INTEGER :: ikt
326      !!--------------------------------------------------------------------
327      !
328      ibld(:,:)   = 0     ; imld(:,:)  = 0
329      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp
330      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp
331      zvw0(:,:)   = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:)      = 0._wp ; zwb0(:,:)      = 0._wp
332      zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:)     = 0._wp ; zwb_ent(:,:)   = 0._wp
333      zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
334      zhol(:,:)   = 0._wp
335      lconv(:,:)  = .FALSE.
336      ! mixed layer
337      ! no initialization of zhbl or zhml (or zdh?)
338      zhbl(:,:)    = 1._wp ; zhml(:,:)    = 1._wp ; zdh(:,:)      = 1._wp ; zdhdt(:,:)   = 0._wp
339      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp
340      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp
341
342      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp
343      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp
344      zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp
345      zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:)  = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp
346      zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp
347      !
348      zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
349      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
350      !
351      zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp
352
353      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed
354         zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp
355         zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp
356         zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp
357         zhmle(:,:) = 0._wp  ; zmld(:,:) = 0._wp
358      ENDIF
359      zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt
360
361      ! Flux-Gradient arrays.
362      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp
363      zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp
364      zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp
365      zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp
366      zhbl_t(:,:)     = 0._wp ; zdhdt(:,:)      = 0._wp
367
368      zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
369      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp
370
371      zdhdt_2(:,:) = 0._wp
372      ! hbl = MAX(hbl,epsln)
373      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
374      ! Calculate boundary layer scales
375      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
376
377      ! Assume two-band radiation model for depth of OSBL
378     zz0 =       rn_abs       ! surface equi-partition in 2-bands
379     zz1 =  1. - rn_abs
380     DO jj = 2, jpjm1
381        DO ji = 2, jpim1
382           ! Surface downward irradiance (so always +ve)
383           zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp
384           ! Downwards irradiance at base of boundary layer
385           zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
386           ! Downwards irradiance averaged over depth of the OSBL
387           zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
388                 &                         + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
389        END DO
390     END DO
391     ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
392     DO jj = 2, jpjm1
393        DO ji = 2, jpim1
394           zthermal = rab_n(ji,jj,1,jp_tem)
395           zbeta    = rab_n(ji,jj,1,jp_sal)
396           ! Upwards surface Temperature flux for non-local term
397           zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
398           ! Upwards surface salinity flux for non-local term
399           zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
400           ! Non radiative upwards surface buoyancy flux
401           zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj)
402           ! turbulent heat flux averaged over depth of OSBL
403           zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
404           ! turbulent salinity flux averaged over depth of the OBSL
405           zwsav(ji,jj) = 0.5 * zws0(ji,jj)
406           ! turbulent buoyancy flux averaged over the depth of the OBSBL
407           zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj)
408           ! Surface upward velocity fluxes
409           zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
410           zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
411           ! Friction velocity (zustar), at T-point : LMD94 eq. 2
412           zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
413           zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
414           zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
415        END DO
416     END DO
417     ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
418     SELECT CASE (nn_osm_wave)
419     ! Assume constant La#=0.3
420     CASE(0)
421        DO jj = 2, jpjm1
422           DO ji = 2, jpim1
423              zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
424              zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
425              zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 )
426              ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
427           END DO
428        END DO
429     ! Assume Pierson-Moskovitz wind-wave spectrum
430     CASE(1)
431        DO jj = 2, jpjm1
432           DO ji = 2, jpim1
433              ! Use wind speed wndm included in sbc_oce module
434              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
435              dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
436           END DO
437        END DO
438     ! Use ECMWF wave fields as output from SBCWAVE
439     CASE(2)
440        zfac =  2.0_wp * rpi / 16.0_wp
441        DO jj = 2, jpjm1
442           DO ji = 2, jpim1
443              ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
444              !    The coefficient 0.8 gives La=0.3  in this situation.
445              ! It could represent the effects of the spread of wave directions
446              ! around the mean wind. The effect of this adjustment needs to be tested.
447              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)
448              zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8)
449              dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
450           END DO
451        END DO
452     END SELECT
453
454     ! Langmuir velocity scale (zwstrl), La # (zla)
455     ! mixed scale (zvstr), convective velocity scale (zwstrc)
456     DO jj = 2, jpjm1
457        DO ji = 2, jpim1
458           ! Langmuir velocity scale (zwstrl), at T-point
459           zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
460           zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)
461           IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))
462           ! Velocity scale that tends to zustar for large Langmuir numbers
463           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + &
464                & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
465
466           ! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
467           ! Note zustke and zwstrl are not amended.
468           !
469           ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
470           IF ( zwbav(ji,jj) > 0.0) THEN
471              zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
472              zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
473              lconv(ji,jj) = .TRUE.
474           ELSE
475              zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln )
476              lconv(ji,jj) = .FALSE.
477           ENDIF
478        END DO
479     END DO
480
481     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
482     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
483     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
484     ! BL must be always 4 levels deep.
485     ! For calculation of lateral buoyancy gradients for FK in
486     ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must
487     ! previously exist for hbl also.
488      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) )
489      ibld(:,:) = 4
490      DO jk = 5, jpkm1
491         DO jj = 1, jpj
492            DO ji = 1, jpi
493               IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
494                  ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
495               ENDIF
496            END DO
497         END DO
498      END DO
499
500      DO jj = 2, jpjm1
501         DO ji = 2, jpim1
502            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
503            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 ))
504            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
505         END DO
506      END DO
507      ! Averages over well-mixed and boundary layer
508      ! Alan: do we need zb_nl?, zb_ml?
509      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl)
510      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
511! External gradient
512      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )
513
514
515      IF ( ln_osm_mle ) THEN
516         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )
517         CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )
518       ENDIF
519
520! Rate of change of hbl
521      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
522      ! Calculate averages over depth of boundary layer
523         DO jj = 2, jpjm1
524            DO ji = 2, jpim1
525               zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it
526               ! adjustment to represent limiting by ocean bottom
527               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:))
528            END DO
529         END DO
530
531      imld(:,:) = ibld(:,:)           ! use imld to hold previous blayer index
532      ibld(:,:) = 4
533
534      DO jk = 4, jpkm1
535         DO jj = 2, jpjm1
536            DO ji = 2, jpim1
537               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
538                  ibld(ji,jj) = jk
539               ENDIF
540            END DO
541         END DO
542      END DO
543
544!
545! Step through model levels taking account of buoyancy change to determine the effect on dhdt
546!
547      CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
548      ! Alan: do we need zb_ml?
549      CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )
550!
551!
552      CALL zdf_osm_pycnocline_thickness( dh, zdh )
553      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms.
554!
555    ! Average over the depth of the mixed layer in the convective boundary layer
556    ! Alan: do we need zb_ml?
557    CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )
558    ! rotate mean currents and changes onto wind align co-ordinates
559    !
560    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
561    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
562     zuw_bse = 0._wp
563     zvw_bse = 0._wp
564     zwth_ent = 0._wp
565     zws_ent = 0._wp
566     DO jj = 2, jpjm1
567        DO ji = 2, jpim1
568           IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN
569              IF ( lconv(ji,jj) ) THEN
570             zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + &
571                      &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ &
572                      &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln)
573            zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ &
574                      &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj))
575                 IF ( zdb_ml(ji,jj) > 0._wp ) THEN
576                    zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
577                    zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
578                 ENDIF
579              ELSE
580                 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )
581                 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )
582              ENDIF
583           ENDIF
584        END DO
585     END DO
586
587      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
588      !  Pycnocline gradients for scalars and velocity
589      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
590
591      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )
592      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc )
593      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )
594       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
595       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
596       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
597
598       DO jj = 2, jpjm1
599          DO ji = 2, jpim1
600             IF ( lconv(ji,jj) ) THEN
601               zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
602               zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
603               zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
604               zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
605               zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third
606               zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln )
607             ELSE
608               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
609               zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
610             END IF
611          END DO
612       END DO
613!
614       DO jj = 2, jpjm1
615          DO ji = 2, jpim1
616             IF ( lconv(ji,jj) ) THEN
617                DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
618                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
619                    !
620                    zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5
621                    !
622                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) &
623                         &            *                                      ( 1.0 -               0.5 * zznd_ml**2 )
624                END DO
625                ! pycnocline - if present linear profile
626                IF ( zdh(ji,jj) > 0._wp ) THEN
627                   zgamma_b = 6.0
628                   DO jk = imld(ji,jj)+1 , ibld(ji,jj)
629                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
630                       !
631                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
632                       !
633                       zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
634                   END DO
635                   IF ( ibld_ext == 0 ) THEN
636                       zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
637                       zviscos(ji,jj,ibld(ji,jj)) = 0._wp
638                   ELSE
639                       zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )
640                       zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )
641                   ENDIF
642                ENDIF
643                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out
644                zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6)
645                zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5)
646                ! could be taken out, take account of entrainment represents as a diffusivity
647                ! should remove w from here, represents entrainment
648             ELSE
649             ! stable conditions
650                DO jk = 2, ibld(ji,jj)
651                   zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
652                   zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
653                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
654                END DO
655
656                IF ( ibld_ext == 0 ) THEN
657                   zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
658                   zviscos(ji,jj,ibld(ji,jj)) = 0._wp
659                ELSE
660                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))
661                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))
662                ENDIF
663             ENDIF   ! end if ( lconv )
664             !
665          END DO  ! end of ji loop
666       END DO     ! end of jj loop
667
668       !
669       ! calculate non-gradient components of the flux-gradient relationships
670       !
671! Stokes term in scalar flux, flux-gradient relationship
672       WHERE ( lconv )
673          zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln)
674          !
675          zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
676       ELSEWHERE
677          zsc_wth_1 = 2.0 * zwthav
678          !
679          zsc_ws_1 = 2.0 * zwsav
680       ENDWHERE
681
682
683       DO jj = 2, jpjm1
684          DO ji = 2, jpim1
685            IF ( lconv(ji,jj) ) THEN
686              DO jk = 2, imld(ji,jj)
687                 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
688                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
689                 !
690                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
691              END DO ! end jk loop
692            ELSE     ! else for if (lconv)
693 ! Stable conditions
694               DO jk = 2, ibld(ji,jj)
695                  zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
696                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
697                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
698                  !
699                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
700                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
701               END DO
702            ENDIF               ! endif for check on lconv
703
704          END DO  ! end of ji loop
705       END DO     ! end of jj loop
706
707! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero)
708       WHERE ( lconv )
709          zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 )
710          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 )
711          zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )
712       ELSEWHERE
713          zsc_uw_1 = zustar**2
714          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln)
715       ENDWHERE
716       IF(ln_dia_osm) THEN
717          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
718          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
719       END IF
720       DO jj = 2, jpjm1
721          DO ji = 2, jpim1
722             IF ( lconv(ji,jj) ) THEN
723                DO jk = 2, imld(ji,jj)
724                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
725                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   &
726                        &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) &
727                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) )
728!
729                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       &
730                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
731                END DO   ! end jk loop
732             ELSE
733! Stable conditions
734                DO jk = 2, ibld(ji,jj) ! corrected to ibld
735                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
736                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       &
737                        &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
738                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
739                END DO   ! end jk loop
740             ENDIF
741          END DO  ! ji loop
742       END DO     ! jj loo
743
744! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
745
746       WHERE ( lconv )
747          zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
748          zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
749       ELSEWHERE
750          zsc_wth_1 = 0._wp
751          zsc_ws_1 = 0._wp
752       ENDWHERE
753
754       DO jj = 2, jpjm1
755          DO ji = 2, jpim1
756             IF (lconv(ji,jj) ) THEN
757                DO jk = 2, imld(ji,jj)
758                   zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
759                   ! calculate turbulent length scale
760                   zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
761                        &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) )
762                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
763                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
764                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
765                   ! non-gradient buoyancy terms
766                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
767                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
768                END DO
769             ELSE
770                DO jk = 2, ibld(ji,jj)
771                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
772                   ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
773                END DO
774             ENDIF
775          END DO   ! ji loop
776       END DO      ! jj loop
777
778       WHERE ( lconv )
779          zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
780          zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
781          zsc_vw_1 = 0._wp
782       ELSEWHERE
783         zsc_uw_1 = 0._wp
784         zsc_vw_1 = 0._wp
785       ENDWHERE
786
787       DO jj = 2, jpjm1
788          DO ji = 2, jpim1
789             IF ( lconv(ji,jj) ) THEN
790                DO jk = 2 , imld(ji,jj)
791                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
792                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     &
793                        &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   &
794                        &                                          * zsc_uw_2(ji,jj)                                    )
795                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
796                END DO  ! jk loop
797             ELSE
798             ! stable conditions
799                DO jk = 2, ibld(ji,jj)
800                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
801                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
802                END DO
803             ENDIF
804          END DO        ! ji loop
805       END DO           ! jj loop
806
807       IF(ln_dia_osm) THEN
808          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
809          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
810       END IF
811! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
812
813       WHERE ( lconv )
814          zsc_wth_1 = zwth0
815          zsc_ws_1 = zws0
816       ELSEWHERE
817          zsc_wth_1 = 2.0 * zwthav
818          zsc_ws_1 = zws0
819       ENDWHERE
820
821       DO jj = 2, jpjm1
822          DO ji = 2, jpim1
823            IF ( lconv(ji,jj) ) THEN
824               DO jk = 2, imld(ji,jj)
825                  zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
826                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                &
827                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
828                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
829                       &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) )
830                  !
831                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  &
832                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
833                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
834                       &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) )
835               END DO
836            ELSE
837               DO jk = 2, ibld(ji,jj)
838                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
839                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
840                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
841               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
842                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
843               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
844               END DO
845            ENDIF
846          ENDDO    ! ji loop
847       END DO      ! jj loop
848
849       WHERE ( lconv )
850          zsc_uw_1 = zustar**2
851          zsc_vw_1 = ff_t * zustke * zhml
852       ELSEWHERE
853          zsc_uw_1 = zustar**2
854          zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
855          zsc_vw_1 = ff_t * zustke * zhbl
856          zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1
857       ENDWHERE
858
859       DO jj = 2, jpjm1
860          DO ji = 2, jpim1
861             IF ( lconv(ji,jj) ) THEN
862               DO jk = 2, imld(ji,jj)
863                  zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
864                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
865                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
866                       & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)
867                  !
868                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
869                       & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
870               END DO
871             ELSE
872               DO jk = 2, ibld(ji,jj)
873                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
874                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
875                  IF ( zznd_d <= 2.0 ) THEN
876                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
877                          &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
878                     !
879                  ELSE
880                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
881                          & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
882                     !
883                  ENDIF
884
885                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
886                       & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj)
887                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
888                       & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
889               END DO
890             ENDIF
891          END DO
892       END DO
893
894       IF(ln_dia_osm) THEN
895          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
896          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
897          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
898          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
899          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
900          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
901       END IF
902!
903! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
904
905      DO jj = 2, jpjm1
906         DO ji = 2, jpim1
907            IF ( lconv(ji,jj) ) THEN
908               DO jk = 2, ibld(ji,jj)
909                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
910                  IF ( znd >= 0.0 ) THEN
911                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
912                     ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
913                  ELSE
914                     ghamu(ji,jj,jk) = 0._wp
915                     ghamv(ji,jj,jk) = 0._wp
916                  ENDIF
917               END DO
918            ELSE
919               DO jk = 2, ibld(ji,jj)
920                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
921                  IF ( znd >= 0.0 ) THEN
922                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
923                     ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
924                  ELSE
925                     ghamu(ji,jj,jk) = 0._wp
926                     ghamv(ji,jj,jk) = 0._wp
927                  ENDIF
928               END DO
929            ENDIF
930         END DO
931      END DO
932
933       IF(ln_dia_osm) THEN
934          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
935          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
936       END IF
937      ! pynocline contributions
938       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small
939       zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln )
940       DO jj = 2, jpjm1
941          DO ji = 2, jpim1
942             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
943                DO jk= 2, ibld(ji,jj)
944                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
945                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
946                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
947                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
948                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk)
949                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
950                END DO
951             END IF
952          END DO
953       END DO
954
955! Entrainment contribution.
956
957       DO jj=2, jpjm1
958          DO ji = 2, jpim1
959            IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN
960               DO jk = 1, imld(ji,jj) - 1
961                  znd=gdepw_n(ji,jj,jk) / zhml(ji,jj)
962                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
963                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
964                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd
965                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd
966               END DO
967               DO jk = imld(ji,jj), ibld(ji,jj)
968                  znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
969                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
970                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
971                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd )
972                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd )
973                END DO
974             ENDIF
975
976             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
977             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
978             ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
979             ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
980          END DO       ! ji loop
981       END DO          ! jj loop
982
983       IF(ln_dia_osm) THEN
984          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
985          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
986          IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse )
987          IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse )
988          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
989          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
990          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
991       END IF
992       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
993       ! Need to put in code for contributions that are applied explicitly to
994       ! the prognostic variables
995       !  1. Entrainment flux
996       !
997       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
998
999
1000
1001       ! rotate non-gradient velocity terms back to model reference frame
1002
1003       DO jj = 2, jpjm1
1004          DO ji = 2, jpim1
1005             DO jk = 2, ibld(ji,jj)
1006                ztemp = ghamu(ji,jj,jk)
1007                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
1008                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
1009             END DO
1010          END DO
1011       END DO
1012
1013       IF(ln_dia_osm) THEN
1014          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1015          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
1016          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
1017       END IF
1018
1019! KPP-style Ri# mixing
1020       IF( ln_kpprimix) THEN
1021          DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
1022             DO jj = 1, jpjm1
1023                DO ji = 1, jpim1   ! vector opt.
1024                   z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
1025                        &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
1026                        &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
1027                   z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
1028                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
1029                        &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
1030                END DO
1031             END DO
1032          END DO
1033      !
1034         DO jk = 2, jpkm1
1035            DO jj = 2, jpjm1
1036               DO ji = 2, jpim1   ! vector opt.
1037                  !                                          ! shear prod. at w-point weightened by mask
1038                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
1039                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
1040                  !                                          ! local Richardson number
1041                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
1042                  zfri =  MIN( zri / rn_riinfty , 1.0_wp )
1043                  zfri  = ( 1.0_wp - zfri * zfri )
1044                  zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk)
1045                END DO
1046             END DO
1047          END DO
1048
1049          DO jj = 2, jpjm1
1050             DO ji = 2, jpim1
1051                DO jk = ibld(ji,jj) + 1, jpkm1
1052                   zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
1053                   zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
1054                END DO
1055             END DO
1056          END DO
1057
1058       END IF ! ln_kpprimix = .true.
1059
1060! KPP-style set diffusivity large if unstable below BL
1061       IF( ln_convmix) THEN
1062          DO jj = 2, jpjm1
1063             DO ji = 2, jpim1
1064                DO jk = ibld(ji,jj) + 1, jpkm1
1065                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
1066                END DO
1067             END DO
1068          END DO
1069       END IF ! ln_convmix = .true.
1070
1071
1072 
1073       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing
1074          DO jj = 2 , jpjm1
1075             DO ji = 2, jpim1
1076                IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL.
1077            ! Calculate MLE flux profiles
1078                   ! DO jk = 1, mld_prof(ji,jj)
1079                   !    znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1080                   !    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + &
1081                   !         & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 )
1082                   !    ghams(ji,jj,jk) = ghams(ji,jj,jk) + &
1083                   !         & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 )
1084                   ! END DO
1085            ! Calculate MLE flux contribution from surface fluxes
1086                   DO jk = 1, ibld(ji,jj)
1087                     znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln)
1088                     ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd )
1089                     ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd )
1090                    END DO
1091                    DO jk = 1, mld_prof(ji,jj)
1092                      znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1093                      ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd )
1094                      ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd )
1095                    END DO
1096            ! Viscosity for MLEs
1097                    DO jk = ibld(ji,jj), mld_prof(ji,jj)
1098                      zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) )
1099                    END DO
1100            ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj)
1101                    jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj))
1102                    jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj))
1103                    DO jk = mld_prof(ji,jj),  jl
1104                       zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * &
1105                            & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / &
1106                            & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln))
1107                    END DO
1108                ENDIF
1109            END DO
1110          END DO
1111       ENDIF
1112
1113       IF(ln_dia_osm) THEN
1114          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1115          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
1116          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
1117       END IF
1118
1119
1120       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
1121       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
1122
1123       ! GN 25/8: need to change tmask --> wmask
1124
1125     DO jk = 2, jpkm1
1126         DO jj = 2, jpjm1
1127             DO ji = 2, jpim1
1128                p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1129                p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
1130             END DO
1131         END DO
1132     END DO
1133      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid  (sign unchanged), needed to caclulate gham[uv] on u and v grids
1134     CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1.,   &
1135      &                  ghamu, 'W', 1. , ghamv, 'W', 1. )
1136       DO jk = 2, jpkm1
1137           DO jj = 2, jpjm1
1138               DO ji = 2, jpim1
1139                  ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
1140                     &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1141
1142                  ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
1143                      &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1144
1145                  ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk)
1146                  ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk)
1147               END DO
1148           END DO
1149        END DO
1150        ! Lateral boundary conditions on final outputs for hbl,  on T-grid (sign unchanged)
1151        CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
1152        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged)
1153        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged)
1154        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1.,   &
1155         &                  ghamu, 'U', -1. , ghamv, 'V', -1. )
1156
1157      IF(ln_dia_osm) THEN
1158         SELECT CASE (nn_osm_wave)
1159         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1).
1160         CASE(0:1)
1161            IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift
1162            IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift
1163            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
1164         ! Stokes drift read in from sbcwave  (=2).
1165         CASE(2)
1166            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift
1167            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift
1168            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period
1169            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height
1170            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
1171            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum
1172            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10
1173            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* &
1174                 & SQRT(ut0sd**2 + vt0sd**2 ) )
1175         END SELECT
1176         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )            ! <Tw_NL>
1177         IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams )            ! <Sw_NL>
1178         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )            ! <uw_NL>
1179         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )            ! <vw_NL>
1180         IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 )            ! <Tw_0>
1181         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0>
1182         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth
1183         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k
1184         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base
1185         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base
1186         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base
1187         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base
1188         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base
1189         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth
1190         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth
1191         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth
1192         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points
1193         IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc )         ! convective velocity scale
1194         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale
1195         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale
1196         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale
1197         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir #
1198         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
1199         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
1200         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine
1201         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine
1202         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine
1203         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine
1204         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine
1205         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux
1206         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux
1207         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux
1208         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux
1209         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML
1210
1211         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )               ! FK layer depth
1212         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )               ! FK target layer depth
1213         IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk )         ! FK b flux
1214         IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b )   ! FK b flux averaged over ML
1215         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k
1216         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )            ! FK dtdx at u-pt
1217         IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy )            ! FK dtdy at v-pt
1218         IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx )            ! FK dtdx at u-pt
1219         IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy )            ! FK dsdy at v-pt
1220         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )            ! FK dbdx at u-pt
1221         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )            ! FK dbdy at v-pt
1222         IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1223         IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1224
1225      END IF
1226
1227CONTAINS
1228
1229
1230   ! Alan: do we need zb?
1231   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv )
1232     !!---------------------------------------------------------------------
1233     !!                   ***  ROUTINE zdf_vertical_average  ***
1234     !!
1235     !! ** Purpose : Determines vertical averages from surface to jnlev.
1236     !!
1237     !! ** Method  : Averages are calculated from the surface to jnlev.
1238     !!              The external level used to calculate differences is ibld+ibld_ext
1239     !!
1240     !!----------------------------------------------------------------------
1241
1242        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over.
1243
1244        ! Alan: do we need zb?
1245        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity
1246        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components
1247        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
1248        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components.
1249
1250        INTEGER :: jk, ji, jj
1251        REAL(wp) :: zthick, zthermal, zbeta
1252
1253
1254        zt   = 0._wp
1255        zs   = 0._wp
1256        zu   = 0._wp
1257        zv   = 0._wp
1258        DO jj = 2, jpjm1                                 !  Vertical slab
1259         DO ji = 2, jpim1
1260            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1261            zbeta    = rab_n(ji,jj,1,jp_sal)
1262               ! average over depth of boundary layer
1263            zthick = epsln
1264            DO jk = 2, jnlev_av(ji,jj)
1265               zthick = zthick + e3t_n(ji,jj,jk)
1266               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
1267               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
1268               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) &
1269                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
1270                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
1271               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) &
1272                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
1273                     &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
1274            END DO
1275            zt(ji,jj) = zt(ji,jj) / zthick
1276            zs(ji,jj) = zs(ji,jj) / zthick
1277            zu(ji,jj) = zu(ji,jj) / zthick
1278            zv(ji,jj) = zv(ji,jj) / zthick
1279            ! Alan: do we need zb?
1280            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem)
1281            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal)
1282            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) &
1283                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) )
1284            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) &
1285                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) )
1286            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
1287         END DO
1288        END DO
1289   END SUBROUTINE zdf_osm_vertical_average
1290
1291   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
1292     !!---------------------------------------------------------------------
1293     !!                   ***  ROUTINE zdf_velocity_rotation  ***
1294     !!
1295     !! ** Purpose : Rotates frame of reference of averaged velocity components.
1296     !!
1297     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w
1298     !!
1299     !!----------------------------------------------------------------------
1300
1301        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle
1302        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current
1303        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline
1304
1305        INTEGER :: ji, jj
1306        REAL(wp) :: ztemp
1307
1308        DO jj = 2, jpjm1
1309           DO ji = 2, jpim1
1310              ztemp = zu(ji,jj)
1311              zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
1312              zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
1313              ztemp = zdu(ji,jj)
1314              zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
1315              zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj)
1316           END DO
1317        END DO
1318    END SUBROUTINE zdf_osm_velocity_rotation
1319
1320    SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz )
1321     !!---------------------------------------------------------------------
1322     !!                   ***  ROUTINE zdf_osm_external_gradients  ***
1323     !!
1324     !! ** Purpose : Calculates the gradients below the OSBL
1325     !!
1326     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient.
1327     !!
1328     !!----------------------------------------------------------------------
1329
1330     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy.
1331
1332     INTEGER :: jj, ji, jkb, jkb1
1333     REAL(wp) :: zthermal, zbeta
1334
1335
1336     DO jj = 2, jpjm1
1337        DO ji = 2, jpim1
1338           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1339              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1340              zbeta    = rab_n(ji,jj,1,jp_sal)
1341              jkb = ibld(ji,jj) + ibld_ext
1342              jkb1 = MIN(jkb + 1, mbkt(ji,jj))
1343              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
1344                   &   / e3t_n(ji,jj,ibld(ji,jj))
1345              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
1346                   &   / e3t_n(ji,jj,ibld(ji,jj))
1347              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
1348           END IF
1349        END DO
1350     END DO
1351    END SUBROUTINE zdf_osm_external_gradients
1352
1353    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz )
1354
1355     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline
1356
1357     INTEGER :: jk, jj, ji
1358     REAL(wp) :: ztgrad, zsgrad, zbgrad
1359     REAL(wp) :: zgamma_b_nd, zgamma_c, znd
1360     REAL(wp) :: zzeta_s=0.3, ztmp
1361
1362     DO jj = 2, jpjm1
1363        DO ji = 2, jpim1
1364           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1365              IF ( lconv(ji,jj) ) THEN  ! convective conditions
1366                    IF ( zdbdz_ext(ji,jj) > 0._wp .AND. &
1367                         & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh &
1368                         &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0
1369                       ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
1370                       ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj)
1371                       zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj)
1372                       zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj)
1373                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
1374                       zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
1375                            & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
1376                       DO jk = 2, ibld(ji,jj)+ibld_ext
1377                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp
1378                          IF ( znd <= zzeta_s ) THEN
1379                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * &
1380                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
1381                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * &
1382                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
1383                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * &
1384                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
1385                          ELSE
1386                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1387                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1388                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1389                          ENDIF
1390                       END DO
1391                    ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing
1392              ELSE
1393                 ! stable conditions
1394                 ! if pycnocline profile only defined when depth steady of increasing.
1395                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady.
1396                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1397                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
1398                          ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
1399                          ztgrad = zdt_bl(ji,jj) * ztmp
1400                          zsgrad = zds_bl(ji,jj) * ztmp
1401                          zbgrad = zdb_bl(ji,jj) * ztmp
1402                          DO jk = 2, ibld(ji,jj)
1403                             znd = gdepw_n(ji,jj,jk) * ztmp
1404                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1405                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1406                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1407                          END DO
1408                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
1409                          ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
1410                          ztgrad = zdt_bl(ji,jj) * ztmp
1411                          zsgrad = zds_bl(ji,jj) * ztmp
1412                          zbgrad = zdb_bl(ji,jj) * ztmp
1413                          DO jk = 2, ibld(ji,jj)
1414                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp
1415                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1416                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1417                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1418                          END DO
1419                       ENDIF ! IF (zhol >=0.5)
1420                    ENDIF    ! IF (zdb_bl> 0.)
1421                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
1422              ENDIF          ! IF (lconv)
1423           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1424        END DO
1425     END DO
1426
1427    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
1428
1429    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
1430      !!---------------------------------------------------------------------
1431      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  ***
1432      !!
1433      !! ** Purpose : Calculates velocity shear in the pycnocline
1434      !!
1435      !! ** Method  :
1436      !!
1437      !!----------------------------------------------------------------------
1438
1439      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
1440
1441      INTEGER :: jk, jj, ji
1442      REAL(wp) :: zugrad, zvgrad, znd
1443      REAL(wp) :: zzeta_v = 0.45
1444      !
1445      DO jj = 2, jpjm1
1446         DO ji = 2, jpim1
1447            !
1448            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1449               IF ( lconv (ji,jj) ) THEN
1450                  ! Unstable conditions
1451                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
1452                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
1453                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
1454                  !Alan is this right?
1455                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
1456                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
1457                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
1458                       &      )/ (zdh(ji,jj)  + epsln )
1459                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
1460                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
1461                     IF ( znd <= 0.0 ) THEN
1462                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
1463                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
1464                     ELSE
1465                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
1466                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
1467                     ENDIF
1468                  END DO
1469               ELSE
1470                  ! stable conditions
1471                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
1472                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
1473                  DO jk = 2, ibld(ji,jj)
1474                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1475                     IF ( znd < 1.0 ) THEN
1476                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
1477                     ELSE
1478                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
1479                     ENDIF
1480                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
1481                  END DO
1482               ENDIF
1483               !
1484            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1485         END DO
1486      END DO
1487    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
1488
1489    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
1490     !!---------------------------------------------------------------------
1491     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
1492     !!
1493     !! ** Purpose : Calculates the rate at which hbl changes.
1494     !!
1495     !! ** Method  :
1496     !!
1497     !!----------------------------------------------------------------------
1498
1499    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl
1500
1501    INTEGER :: jj, ji
1502    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert
1503    REAL(wp) :: zvel_max, zwb_min
1504    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
1505    REAL(wp) :: zzeta_m = 0.3
1506    REAL(wp) :: zgamma_c = 2.0
1507    REAL(wp) :: zdhoh = 0.1
1508    REAL(wp) :: alpha_bc = 0.5
1509
1510    DO jj = 2, jpjm1
1511       DO ji = 2, jpim1
1512          IF ( lconv(ji,jj) ) THEN    ! Convective
1513             ! Alan is this right?  Yes, it's a bit different from the previous relationship
1514             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) &
1515             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1516             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
1517             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
1518             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
1519             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
1520             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
1521                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
1522
1523             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
1524                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
1525                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
1526                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1527             !
1528             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj)
1529
1530             IF ( ln_osm_mle ) THEN
1531                  !  zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + &
1532                ! &            ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 )          ! Fox-Kemper buoyancy flux average over OSBL
1533                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
1534                   zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  &
1535                        (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
1536                ELSE
1537                   zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
1538                ENDIF
1539                zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1540                IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
1541                   ! OSBL is deepening, entrainment > restratification
1542                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1543                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1544                   ELSE
1545                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15)
1546                   ENDIF
1547                ELSE
1548                   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
1549                   zdhdt(ji,jj) = - zvel_mle(ji,jj)
1550
1551
1552                ENDIF
1553
1554             ELSE
1555                ! Fox-Kemper not used.
1556
1557                  zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
1558                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1559                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1560                ! added ajgn 23 July as temporay fix
1561
1562             ENDIF
1563
1564             zdhdt_2(ji,jj) = 0._wp
1565
1566                ! commented out ajgn 23 July as temporay fix
1567!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1568! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary.
1569!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1570!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj)
1571!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1572!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj)
1573!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min )
1574!                     ! Alan no idea what this should be?
1575!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) &
1576!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) &
1577!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj)
1578!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
1579!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN
1580!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj)
1581!                 ENDIF
1582            ELSE                        ! Stable
1583                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
1584                zdhdt_2(ji,jj) = 0._wp
1585                IF ( zdhdt(ji,jj) < 0._wp ) THEN
1586                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
1587                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
1588                ELSE
1589                    zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
1590                ENDIF
1591                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert
1592            ENDIF
1593         END DO
1594      END DO
1595    END SUBROUTINE zdf_osm_calculate_dhdt
1596
1597    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
1598     !!---------------------------------------------------------------------
1599     !!                   ***  ROUTINE zdf_osm_timestep_hbl  ***
1600     !!
1601     !! ** Purpose : Increments hbl.
1602     !!
1603     !! ** Method  : If thechange in hbl exceeds one model level the change is
1604     !!              is calculated by moving down the grid, changing the buoyancy
1605     !!              jump. This is to ensure that the change in hbl does not
1606     !!              overshoot a stable layer.
1607     !!
1608     !!----------------------------------------------------------------------
1609
1610
1611    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl.
1612
1613    INTEGER :: jk, jj, ji, jm
1614    REAL(wp) :: zhbl_s, zvel_max, zdb
1615    REAL(wp) :: zthermal, zbeta
1616
1617     DO jj = 2, jpjm1
1618         DO ji = 2, jpim1
1619           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
1620!
1621! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
1622!
1623              zhbl_s = hbl(ji,jj)
1624              jm = imld(ji,jj)
1625              zthermal = rab_n(ji,jj,1,jp_tem)
1626              zbeta = rab_n(ji,jj,1,jp_sal)
1627
1628
1629              IF ( lconv(ji,jj) ) THEN
1630!unstable
1631
1632                 IF( ln_osm_mle ) THEN
1633                    zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1634                 ELSE
1635
1636                    zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
1637                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1638
1639                 ENDIF
1640
1641                 DO jk = imld(ji,jj), ibld(ji,jj)
1642                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
1643                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
1644                         &  0.0 ) + zvel_max
1645
1646
1647                    IF ( ln_osm_mle ) THEN
1648                       zhbl_s = zhbl_s + MIN( &
1649                        & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
1650                        & e3w_n(ji,jj,jm) )
1651                    ELSE
1652                      zhbl_s = zhbl_s + MIN( &
1653                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
1654                        & e3w_n(ji,jj,jm) )
1655                    ENDIF
1656
1657                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1658
1659                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
1660                 END DO
1661                 hbl(ji,jj) = zhbl_s
1662                 ibld(ji,jj) = jm
1663             ELSE
1664! stable
1665                 DO jk = imld(ji,jj), ibld(ji,jj)
1666                    zdb = MAX( &
1667                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&
1668                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
1669                         & 0.0 ) + &
1670             &       2.0 * zvstr(ji,jj)**2 / zhbl_s
1671
1672                    ! Alan is thuis right? I have simply changed hbli to hbl
1673                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
1674                    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) ) ) * &
1675               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
1676                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
1677                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
1678
1679                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1680                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
1681                 END DO
1682             ENDIF   ! IF ( lconv )
1683             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
1684             ibld(ji,jj) = MAX(jm, 4 )
1685           ELSE
1686! change zero or one model level.
1687             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
1688           ENDIF
1689           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
1690         END DO
1691      END DO
1692
1693    END SUBROUTINE zdf_osm_timestep_hbl
1694
1695    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
1696      !!---------------------------------------------------------------------
1697      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  ***
1698      !!
1699      !! ** Purpose : Calculates thickness of the pycnocline
1700      !!
1701      !! ** Method  : The thickness is calculated from a prognostic equation
1702      !!              that relaxes the pycnocine thickness to a diagnostic
1703      !!              value. The time change is calculated assuming the
1704      !!              thickness relaxes exponentially. This is done to deal
1705      !!              with large timesteps.
1706      !!
1707      !!----------------------------------------------------------------------
1708
1709      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness.
1710      !
1711      INTEGER :: jj, ji
1712      INTEGER :: inhml
1713      REAL(wp) :: zari, ztau, zddhdt
1714
1715
1716      DO jj = 2, jpjm1
1717         DO ji = 2, jpim1
1718            IF ( lconv(ji,jj) ) THEN
1719
1720               IF( ln_osm_mle ) THEN
1721                  IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN
1722                     ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
1723                     IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1724                        IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
1725                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1726                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1727                        ELSE                                                     ! unstable
1728                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1729                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1730                        ENDIF
1731                        ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
1732                        zddhdt = zari * hbl(ji,jj)
1733                     ELSE
1734                        ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1735                        zddhdt = 0.2 * hbl(ji,jj)
1736                     ENDIF
1737                  ELSE
1738                     ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1739                     zddhdt = 0.2 * hbl(ji,jj)
1740                  ENDIF
1741               ELSE ! ln_osm_mle
1742                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1743                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
1744                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1745                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1746                     ELSE                                                     ! unstable
1747                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1748                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1749                     ENDIF
1750                     ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
1751                     zddhdt = zari * hbl(ji,jj)
1752                  ELSE
1753                     ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1754                     zddhdt = 0.2 * hbl(ji,jj)
1755                  ENDIF
1756
1757               END IF  ! ln_osm_mle
1758
1759               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1760               ! Alan: this hml is never defined or used
1761            ELSE   ! IF (lconv)
1762               ztau = hbl(ji,jj) / zvstr(ji,jj)
1763               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
1764                  ! boundary layer deepening
1765                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1766                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
1767                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
1768                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
1769                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj)
1770                  ELSE
1771                     zddhdt = 0.2 * hbl(ji,jj)
1772                  ENDIF
1773               ELSE     ! IF(dhdt < 0)
1774                  zddhdt = 0.2 * hbl(ji,jj)
1775               ENDIF    ! IF (dhdt >= 0)
1776               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1777               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt       ! can be a problem with dh>hbl for rapid collapse
1778               ! Alan: this hml is never defined or used -- do we need it?
1779            ENDIF       ! IF (lconv)
1780
1781            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
1782            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
1783            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
1784            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
1785            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1786         END DO
1787      END DO
1788
1789    END SUBROUTINE zdf_osm_pycnocline_thickness
1790
1791
1792   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )
1793      !!----------------------------------------------------------------------
1794      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  ***
1795      !!
1796      !! ** Purpose :   Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.
1797      !!
1798      !! ** Method  :
1799      !!
1800      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
1801      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
1802
1803
1804      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points
1805      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == !
1806      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy
1807
1808      INTEGER  ::   ji, jj, jk          ! dummy loop indices
1809      INTEGER  ::   ii, ij, ik, ikmax   ! local integers
1810      REAL(wp)                         :: zc
1811      REAL(wp)                         :: zN2_c           ! local buoyancy difference from 10m value
1812      REAL(wp), DIMENSION(jpi,jpj)     :: ztm, zsm, zLf_NH, zLf_MH
1813      REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv
1814      REAL(wp), DIMENSION(jpi,jpj)     :: zmld_midu, zmld_midv
1815!!----------------------------------------------------------------------
1816      !
1817      !                                      !==  MLD used for MLE  ==!
1818
1819      mld_prof(:,:)  = nlb10               ! Initialization to the number of w ocean point
1820      zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
1821      zN2_c = grav * rn_osm_mle_rho_c * r1_rau0   ! convert density criteria into N^2 criteria
1822      DO jk = nlb10, jpkm1
1823         DO jj = 1, jpj                ! Mixed layer level: w-level
1824            DO ji = 1, jpi
1825               ikt = mbkt(ji,jj)
1826               zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
1827               IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
1828            END DO
1829         END DO
1830      END DO
1831      DO jj = 1, jpj
1832         DO ji = 1, jpi
1833            mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
1834            zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
1835         END DO
1836      END DO
1837      ! ensure mld_prof .ge. ibld
1838      !
1839      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation
1840      !
1841      ztm(:,:) = 0._wp
1842      zsm(:,:) = 0._wp
1843      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
1844         DO jj = 1, jpj
1845            DO ji = 1, jpi
1846               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
1847               ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem)
1848               zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal)
1849            END DO
1850         END DO
1851      END DO
1852      ! average temperature and salinity.
1853      ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
1854      zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
1855      ! calculate horizontal gradients at u & v points
1856
1857      DO jj = 2, jpjm1
1858         DO ji = 1, jpim1
1859            zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
1860            zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
1861            zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))
1862            ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )
1863            ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )
1864         END DO
1865      END DO
1866
1867      DO jj = 1, jpjm1
1868         DO ji = 2, jpim1
1869            zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
1870            zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
1871            zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))
1872            ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )
1873            ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )
1874         END DO
1875      END DO
1876
1877      CALL eos_rab(ztsm_midu, zmld_midu, zabu)
1878      CALL eos_rab(ztsm_midv, zmld_midv, zabv)
1879
1880      DO jj = 2, jpjm1
1881         DO ji = 1, jpim1
1882            dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
1883         END DO
1884      END DO
1885      DO jj = 1, jpjm1
1886         DO ji = 2, jpim1
1887            dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
1888         END DO
1889      END DO
1890
1891 END SUBROUTINE zdf_osm_zmld_horizontal_gradients
1892  SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )
1893      !!----------------------------------------------------------------------
1894      !!                  ***  ROUTINE zdf_osm_mle_parameters  ***
1895      !!
1896      !! ** Purpose :   Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.
1897      !!
1898      !! ** Method  :
1899      !!
1900      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
1901      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
1902
1903      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zwb_fk, zvel_mle, zdiff_mle
1904      INTEGER  ::   ji, jj, jk          ! dummy loop indices
1905      INTEGER  ::   ii, ij, ik  ! local integers
1906      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
1907      REAL(wp) ::   zdb_mle, ztmp, zdbds_mle
1908
1909      mld_prof(:,:) = 4
1910      DO jk = 5, jpkm1
1911        DO jj = 2, jpjm1
1912          DO ji = 2, jpim1
1913            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
1914          END DO
1915        END DO
1916      END DO
1917      ! DO jj = 2, jpjm1
1918      !    DO ji = 1, jpim1
1919      !       zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
1920      !    END DO
1921      !  END DO
1922   ! Timestep mixed layer eddy depth.
1923      DO jj = 2, jpjm1
1924        DO ji = 2, jpim1
1925          zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG
1926          IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN
1927             hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh )  ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy.
1928          ELSE
1929            ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10
1930             ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN
1931             !   hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau
1932             ! ELSE
1933             !   hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD
1934             ! ENDIF
1935             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
1936               hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau
1937             ELSE
1938               hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD
1939             ENDIF
1940          ENDIF
1941          hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj))
1942          hmle(ji,jj) = MIN(hmle(ji,jj), 1.2*zmld(ji,jj)) 
1943        END DO
1944      END DO
1945
1946      mld_prof = 4
1947      DO jk = 5, jpkm1
1948        DO jj = 2, jpjm1
1949          DO ji = 2, jpim1
1950            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
1951          END DO
1952        END DO
1953      END DO
1954      DO jj = 2, jpjm1
1955         DO ji = 2, jpim1
1956            zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj))
1957         END DO
1958       END DO
1959   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
1960
1961      DO jj = 2, jpjm1
1962        DO ji = 2, jpim1
1963          IF ( lconv(ji,jj) ) THEN
1964             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
1965             ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) &
1966             !      & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) )
1967             ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) &
1968             !      & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) )
1969             zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &
1970                  & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
1971             zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle
1972      ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
1973             zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
1974             zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf
1975          ENDIF
1976        END DO
1977      END DO
1978END SUBROUTINE zdf_osm_mle_parameters
1979
1980END SUBROUTINE zdf_osm
1981
1982
1983   SUBROUTINE zdf_osm_init
1984     !!----------------------------------------------------------------------
1985     !!                  ***  ROUTINE zdf_osm_init  ***
1986     !!
1987     !! ** Purpose :   Initialization of the vertical eddy diffivity and
1988     !!      viscosity when using a osm turbulent closure scheme
1989     !!
1990     !! ** Method  :   Read the namosm namelist and check the parameters
1991     !!      called at the first timestep (nit000)
1992     !!
1993     !! ** input   :   Namlist namosm
1994     !!----------------------------------------------------------------------
1995     INTEGER  ::   ios            ! local integer
1996     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1997     REAL z1_t2
1998     !!
1999     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
2000          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
2001          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
2002          & ,ln_osm_mle
2003! Namelist for Fox-Kemper parametrization.
2004      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
2005           & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau
2006
2007     !!----------------------------------------------------------------------
2008     !
2009     REWIND( numnam_ref )              ! Namelist namzdf_osm in reference namelist : Osmosis ML model
2010     READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
2011901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
2012
2013     REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
2014     READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
2015902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
2016     IF(lwm) WRITE ( numond, namzdf_osm )
2017
2018     IF(lwp) THEN                    ! Control print
2019        WRITE(numout,*)
2020        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
2021        WRITE(numout,*) '~~~~~~~~~~~~'
2022        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters'
2023        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la
2024        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle
2025        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
2026        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
2027        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes
2028        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
2029        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave
2030        SELECT CASE (nn_osm_wave)
2031        CASE(0)
2032           WRITE(numout,*) '     calculated assuming constant La#=0.3'
2033        CASE(1)
2034           WRITE(numout,*) '     calculated from Pierson Moskowitz wind-waves'
2035        CASE(2)
2036           WRITE(numout,*) '     calculated from ECMWF wave fields'
2037        END SELECT
2038        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
2039        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix
2040        WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
2041        WRITE(numout,*) '     maximum shear diffusivity at Rig = 0    (m2/s) rn_difri = ', rn_difri
2042        WRITE(numout,*) '     Use large mixing below BL when unstable       ln_convmix = ', ln_convmix
2043        WRITE(numout,*) '     diffusivity when unstable below BL     (m2/s) rn_difconv = ', rn_difconv
2044     ENDIF
2045
2046     !                              ! allocate zdfosm arrays
2047     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
2048
2049
2050     IF( ln_osm_mle ) THEN
2051! Initialise Fox-Kemper parametrization
2052         REWIND( numnam_ref )              ! Namelist namosm_mle in reference namelist : Tracer advection scheme
2053         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
2054903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist')
2055
2056         REWIND( numnam_cfg )              ! Namelist namosm_mle in configuration namelist : Tracer advection scheme
2057         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
2058904      IF( ios >  0 )   CALL ctl_nam ( ios , 'namosm_mle in configuration namelist')
2059         IF(lwm) WRITE ( numond, namosm_mle )
2060
2061         IF(lwp) THEN                     ! Namelist print
2062            WRITE(numout,*)
2063            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
2064            WRITE(numout,*) '~~~~~~~~~~~~~'
2065            WRITE(numout,*) '   Namelist namosm_mle : '
2066            WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_osm_mle    = ', nn_osm_mle
2067            WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_osm_mle_ce    = ', rn_osm_mle_ce
2068            WRITE(numout,*) '         scale of ML front (ML radius of deformation) (nn_osm_mle=0)      rn_osm_mle_lf     = ', rn_osm_mle_lf, 'm'
2069            WRITE(numout,*) '         maximum time scale of MLE                    (nn_osm_mle=0)      rn_osm_mle_time   = ', rn_osm_mle_time, 's'
2070            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg'
2071            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c
2072            WRITE(numout,*) '         Threshold used to define ML for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s'
2073            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's'
2074         ENDIF         !
2075     ENDIF
2076      !
2077      IF(lwp) THEN
2078         WRITE(numout,*)
2079         IF( ln_osm_mle ) THEN
2080            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
2081            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
2082            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation'
2083         ELSE
2084            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
2085         ENDIF
2086      ENDIF
2087      !
2088      IF( ln_osm_mle ) THEN                ! MLE initialisation
2089         !
2090         rb_c = grav * rn_osm_mle_rho_c /rau0        ! Mixed Layer buoyancy criteria
2091         IF(lwp) WRITE(numout,*)
2092         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
2093         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
2094         !
2095         IF( nn_osm_mle == 0 ) THEN           ! MLE array allocation & initialisation            !
2096!
2097         ELSEIF( nn_osm_mle == 1 ) THEN           ! MLE array allocation & initialisation
2098            rc_f = rn_osm_mle_ce/ (  5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat )  )
2099            !
2100         ENDIF
2101         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
2102         z1_t2 = 2.e-5
2103         do jj=1,jpj
2104            do ji = 1,jpi
2105               r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
2106            end do
2107         end do
2108         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
2109         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
2110         !
2111      ENDIF
2112
2113     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle
2114
2115
2116     IF( ln_zdfddm) THEN
2117        IF(lwp) THEN
2118           WRITE(numout,*)
2119           WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
2120           WRITE(numout,*) '    CAUTION : done in routine zdfosm, not in routine zdfddm '
2121        ENDIF
2122     ENDIF
2123
2124
2125     !set constants not in namelist
2126     !-----------------------------
2127
2128     IF(lwp) THEN
2129        WRITE(numout,*)
2130     ENDIF
2131
2132     IF (nn_osm_wave == 0) THEN
2133        dstokes(:,:) = rn_osm_dstokes
2134     END IF
2135
2136     ! Horizontal average : initialization of weighting arrays
2137     ! -------------------
2138
2139     SELECT CASE ( nn_ave )
2140
2141     CASE ( 0 )                ! no horizontal average
2142        IF(lwp) WRITE(numout,*) '          no horizontal average on avt'
2143        IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
2144        ! weighting mean arrays etmean
2145        !           ( 1  1 )
2146        ! avt = 1/4 ( 1  1 )
2147        !
2148        etmean(:,:,:) = 0.e0
2149
2150        DO jk = 1, jpkm1
2151           DO jj = 2, jpjm1
2152              DO ji = 2, jpim1   ! vector opt.
2153                 etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
2154                      &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
2155                      &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
2156              END DO
2157           END DO
2158        END DO
2159
2160     CASE ( 1 )                ! horizontal average
2161        IF(lwp) WRITE(numout,*) '          horizontal average on avt'
2162        ! weighting mean arrays etmean
2163        !           ( 1/2  1  1/2 )
2164        ! avt = 1/8 ( 1    2  1   )
2165        !           ( 1/2  1  1/2 )
2166        etmean(:,:,:) = 0.e0
2167
2168        DO jk = 1, jpkm1
2169           DO jj = 2, jpjm1
2170              DO ji = 2, jpim1   ! vector opt.
2171                 etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
2172                      & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
2173                      &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
2174                      &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
2175                      &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
2176                      &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
2177              END DO
2178           END DO
2179        END DO
2180
2181     CASE DEFAULT
2182        WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
2183        CALL ctl_stop( ctmp1 )
2184
2185     END SELECT
2186
2187     ! Initialization of vertical eddy coef. to the background value
2188     ! -------------------------------------------------------------
2189     DO jk = 1, jpk
2190        avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
2191     END DO
2192
2193     ! zero the surface flux for non local term and osm mixed layer depth
2194     ! ------------------------------------------------------------------
2195     ghamt(:,:,:) = 0.
2196     ghams(:,:,:) = 0.
2197     ghamu(:,:,:) = 0.
2198     ghamv(:,:,:) = 0.
2199     !
2200     IF( lwxios ) THEN
2201        CALL iom_set_rstw_var_active('wn')
2202        CALL iom_set_rstw_var_active('hbl')
2203        CALL iom_set_rstw_var_active('dh')
2204        IF( ln_osm_mle ) THEN
2205            CALL iom_set_rstw_var_active('hmle')
2206        END IF
2207     ENDIF
2208   END SUBROUTINE zdf_osm_init
2209
2210
2211   SUBROUTINE osm_rst( kt, cdrw )
2212     !!---------------------------------------------------------------------
2213     !!                   ***  ROUTINE osm_rst  ***
2214     !!
2215     !! ** Purpose :   Read or write BL fields in restart file
2216     !!
2217     !! ** Method  :   use of IOM library. If the restart does not contain
2218     !!                required fields, they are recomputed from stratification
2219     !!----------------------------------------------------------------------
2220
2221     INTEGER, INTENT(in) :: kt
2222     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
2223
2224     INTEGER ::   id1, id2, id3   ! iom enquiry index
2225     INTEGER  ::   ji, jj, jk     ! dummy loop indices
2226     INTEGER  ::   iiki, ikt ! local integer
2227     REAL(wp) ::   zhbf           ! tempory scalars
2228     REAL(wp) ::   zN2_c           ! local scalar
2229     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
2230     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
2231     !!----------------------------------------------------------------------
2232     !
2233     !!-----------------------------------------------------------------------------
2234     ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
2235     !!-----------------------------------------------------------------------------
2236     IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
2237        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. )
2238        IF( id1 > 0 ) THEN                       ! 'wn' exists; read
2239           CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios )
2240           WRITE(numout,*) ' ===>>>> :  wn read from restart file'
2241        ELSE
2242           wn(:,:,:) = 0._wp
2243           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
2244        END IF
2245
2246        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
2247        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. )
2248        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return
2249           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
2250           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  )
2251           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file'
2252           IF( ln_osm_mle ) THEN
2253              id3 = iom_varid( numror, 'hmle'   , ldstop = .FALSE. )
2254              IF( id3 > 0) THEN
2255                 CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios )
2256                 WRITE(numout,*) ' ===>>>> :  hmle read from restart file'
2257              ELSE
2258                 WRITE(numout,*) ' ===>>>> :  hmle not found, set to hbl'
2259                 hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
2260              END IF
2261           END IF
2262           RETURN
2263        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate
2264           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
2265        END IF
2266     END IF
2267
2268     !!-----------------------------------------------------------------------------
2269     ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
2270     !!-----------------------------------------------------------------------------
2271     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
2272        IF(lwp) WRITE(numout,*) '---- osm-rst ----'
2273         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn,   ldxios = lwxios )
2274         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl,  ldxios = lwxios )
2275         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh,   ldxios = lwxios )
2276         IF( ln_osm_mle ) THEN
2277            CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios )
2278         END IF
2279        RETURN
2280     END IF
2281
2282     !!-----------------------------------------------------------------------------
2283     ! Getting hbl, no restart file with hbl, so calculate from surface stratification
2284     !!-----------------------------------------------------------------------------
2285     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
2286     ! w-level of the mixing and mixed layers
2287     CALL eos_rab( tsn, rab_n )
2288     CALL bn2(tsn, rab_n, rn2)
2289     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
2290     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
2291     zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
2292     !
2293     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
2294     DO jk = 1, jpkm1
2295        DO jj = 1, jpj              ! Mixed layer level: w-level
2296           DO ji = 1, jpi
2297              ikt = mbkt(ji,jj)
2298              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
2299              IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
2300           END DO
2301        END DO
2302     END DO
2303     !
2304     DO jj = 1, jpj
2305        DO ji = 1, jpi
2306           iiki = MAX(4,imld_rst(ji,jj))
2307           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth
2308           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth
2309        END DO
2310     END DO
2311
2312     WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
2313
2314     IF( ln_osm_mle ) THEN
2315        hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
2316        WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
2317     END IF
2318
2319     wn(:,:,:) = 0._wp
2320     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
2321   END SUBROUTINE osm_rst
2322
2323
2324   SUBROUTINE tra_osm( kt )
2325      !!----------------------------------------------------------------------
2326      !!                  ***  ROUTINE tra_osm  ***
2327      !!
2328      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
2329      !!
2330      !! ** Method  :   ???
2331      !!----------------------------------------------------------------------
2332      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
2333      !!----------------------------------------------------------------------
2334      INTEGER, INTENT(in) :: kt
2335      INTEGER :: ji, jj, jk
2336      !
2337      IF( kt == nit000 ) THEN
2338         IF(lwp) WRITE(numout,*)
2339         IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
2340         IF(lwp) WRITE(numout,*) '~~~~~~~   '
2341      ENDIF
2342
2343      IF( l_trdtra )   THEN                    !* Save ta and sa trends
2344         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
2345         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
2346      ENDIF
2347
2348      DO jk = 1, jpkm1
2349         DO jj = 2, jpjm1
2350            DO ji = 2, jpim1
2351               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
2352                  &                 - (  ghamt(ji,jj,jk  )  &
2353                  &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
2354               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
2355                  &                 - (  ghams(ji,jj,jk  )  &
2356                  &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
2357            END DO
2358         END DO
2359      END DO
2360
2361      ! save the non-local tracer flux trends for diagnostics
2362      IF( l_trdtra )   THEN
2363         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
2364         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
2365
2366         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt )
2367         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds )
2368         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
2369      ENDIF
2370
2371      IF(ln_ctl) THEN
2372         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   &
2373         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
2374      ENDIF
2375      !
2376   END SUBROUTINE tra_osm
2377
2378
2379   SUBROUTINE trc_osm( kt )          ! Dummy routine
2380      !!----------------------------------------------------------------------
2381      !!                  ***  ROUTINE trc_osm  ***
2382      !!
2383      !! ** Purpose :   compute and add to the passive tracer trend the non-local
2384      !!                 passive tracer flux
2385      !!
2386      !!
2387      !! ** Method  :   ???
2388      !!----------------------------------------------------------------------
2389      !
2390      !!----------------------------------------------------------------------
2391      INTEGER, INTENT(in) :: kt
2392      WRITE(*,*) 'trc_osm: Not written yet', kt
2393   END SUBROUTINE trc_osm
2394
2395
2396   SUBROUTINE dyn_osm( kt )
2397      !!----------------------------------------------------------------------
2398      !!                  ***  ROUTINE dyn_osm  ***
2399      !!
2400      !! ** Purpose :   compute and add to the velocity trend the non-local flux
2401      !! copied/modified from tra_osm
2402      !!
2403      !! ** Method  :   ???
2404      !!----------------------------------------------------------------------
2405      INTEGER, INTENT(in) ::   kt   !
2406      !
2407      INTEGER :: ji, jj, jk   ! dummy loop indices
2408      !!----------------------------------------------------------------------
2409      !
2410      IF( kt == nit000 ) THEN
2411         IF(lwp) WRITE(numout,*)
2412         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
2413         IF(lwp) WRITE(numout,*) '~~~~~~~   '
2414      ENDIF
2415      !code saving tracer trends removed, replace with trdmxl_oce
2416
2417      DO jk = 1, jpkm1           ! add non-local u and v fluxes
2418         DO jj = 2, jpjm1
2419            DO ji = 2, jpim1
2420               ua(ji,jj,jk) =  ua(ji,jj,jk)                      &
2421                  &                 - (  ghamu(ji,jj,jk  )  &
2422                  &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
2423               va(ji,jj,jk) =  va(ji,jj,jk)                      &
2424                  &                 - (  ghamv(ji,jj,jk  )  &
2425                  &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
2426            END DO
2427         END DO
2428      END DO
2429      !
2430      ! code for saving tracer trends removed
2431      !
2432   END SUBROUTINE dyn_osm
2433
2434   !!======================================================================
2435
2436END MODULE zdfosm
Note: See TracBrowser for help on using the repository browser.