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_FKOSM/src/OCE/ZDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfosm.F90 @ 12322

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

Update with George's latest changes for restartability and reproducibility

File size: 127.1 KB
Line 
1MODULE zdfosm
2   !!======================================================================
3   !!                       ***  MODULE  zdfosm  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the OSMOSIS
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !!  History : NEMO 4.0  ! A. Grant, G. Nurser
8   !! 15/03/2017  Changed calculation of pycnocline thickness in unstable conditions and stable conditions AG
9   !! 15/03/2017  Calculation of pycnocline gradients for stable conditions changed. Pycnocline gradients now depend on stability of the OSBL. A.G
10   !! 06/06/2017  (1) Checks on sign of buoyancy jump in calculation of  OSBL depth.  A.G.
11   !!             (2) Removed variable zbrad0, zbradh and zbradav since they are not used.
12   !!             (3) Approximate treatment for shear turbulence.
13   !!                        Minimum values for zustar and zustke.
14   !!                        Add velocity scale, zvstr, that tends to zustar for large Langmuir numbers.
15   !!                        Limit maximum value for Langmuir number.
16   !!                        Use zvstr in definition of stability parameter zhol.
17   !!             (4) Modified parametrization of entrainment flux, changing original coefficient 0.0485 for Langmuir contribution to 0.135 * zla
18   !!             (5) For stable boundary layer add factor that depends on length of timestep to 'slow' collapse and growth. Make sure buoyancy jump not negative.
19   !!             (6) For unstable conditions when growth is over multiple levels, limit change to maximum of one level per cycle through loop.
20   !!             (7) Change lower limits for loops that calculate OSBL averages from 1 to 2. Large gradients between levels 1 and 2 can cause problems.
21   !!             (8) Change upper limits from ibld-1 to ibld.
22   !!             (9) Calculation of pycnocline thickness in unstable conditions. Check added to ensure that buoyancy jump is positive before calculating Ri.
23   !!            (10) Thickness of interface layer at base of the stable OSBL set by Richardson number. Gives continuity in transition from unstable OSBL.
24   !!            (11) Checks that buoyancy jump is poitive when calculating pycnocline profiles.
25   !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity.
26   !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information
27   !!            (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence.
28   !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added.
29   !!            (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out)
30   !!            (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out)
31   !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made,
32   !!                  (a) Pycnocline temperature and salinity profies changed for unstable layers
33   !!                  (b) The stable OSBL depth parametrization changed.
34   !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code.
35   !! 23/05/19   (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1
36   !!----------------------------------------------------------------------
37
38   !!----------------------------------------------------------------------
39   !!   'ln_zdfosm'                                             OSMOSIS scheme
40   !!----------------------------------------------------------------------
41   !!   zdf_osm       : update momentum and tracer Kz from osm scheme
42   !!   zdf_osm_init  : initialization, namelist read, and parameters control
43   !!   osm_rst       : read (or initialize) and write osmosis restart fields
44   !!   tra_osm       : compute and add to the T & S trend the non-local flux
45   !!   trc_osm       : compute and add to the passive tracer trend the non-local flux (TBD)
46   !!   dyn_osm       : compute and add to u & v trensd the non-local flux
47   !!
48   !! Subroutines in revised code.
49   !!----------------------------------------------------------------------
50   USE oce            ! ocean dynamics and active tracers
51                      ! uses wn from previous time step (which is now wb) to calculate hbl
52   USE dom_oce        ! ocean space and time domain
53   USE zdf_oce        ! ocean vertical physics
54   USE sbc_oce        ! surface boundary condition: ocean
55   USE sbcwave        ! surface wave parameters
56   USE phycst         ! physical constants
57   USE eosbn2         ! equation of state
58   USE traqsr         ! details of solar radiation absorption
59   USE zdfddm         ! double diffusion mixing (avs array)
60   USE 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
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                       ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj)
1370                       zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj)
1371                       zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj)
1372                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1373                       zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
1374                            & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
1375                       DO jk = 2, ibld(ji,jj)+ibld_ext
1376                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1377                          IF ( znd <= zzeta_s ) THEN
1378                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * &
1379                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
1380                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * &
1381                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
1382                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * &
1383                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
1384                          ELSE
1385                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1386                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1387                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1388                          ENDIF
1389                       END DO
1390                    ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing
1391              ELSE
1392                 ! stable conditions
1393                 ! if pycnocline profile only defined when depth steady of increasing.
1394                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady.
1395                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1396                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
1397                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
1398                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
1399                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
1400                          DO jk = 2, ibld(ji,jj)
1401                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1402                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1403                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1404                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1405                          END DO
1406                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
1407                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
1408                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
1409                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
1410                          DO jk = 2, ibld(ji,jj)
1411                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
1412                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1413                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1414                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1415                          END DO
1416                       ENDIF ! IF (zhol >=0.5)
1417                    ENDIF    ! IF (zdb_bl> 0.)
1418                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
1419              ENDIF          ! IF (lconv)
1420           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1421        END DO
1422     END DO
1423
1424    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
1425
1426    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
1427      !!---------------------------------------------------------------------
1428      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  ***
1429      !!
1430      !! ** Purpose : Calculates velocity shear in the pycnocline
1431      !!
1432      !! ** Method  :
1433      !!
1434      !!----------------------------------------------------------------------
1435
1436      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
1437
1438      INTEGER :: jk, jj, ji
1439      REAL(wp) :: zugrad, zvgrad, znd
1440      REAL(wp) :: zzeta_v = 0.45
1441      !
1442      DO jj = 2, jpjm1
1443         DO ji = 2, jpim1
1444            !
1445            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1446               IF ( lconv (ji,jj) ) THEN
1447                  ! Unstable conditions
1448                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
1449                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
1450                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
1451                  !Alan is this right?
1452                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
1453                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
1454                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
1455                       &      )/ (zdh(ji,jj)  + epsln )
1456                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
1457                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
1458                     IF ( znd <= 0.0 ) THEN
1459                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
1460                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
1461                     ELSE
1462                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
1463                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
1464                     ENDIF
1465                  END DO
1466               ELSE
1467                  ! stable conditions
1468                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
1469                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
1470                  DO jk = 2, ibld(ji,jj)
1471                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1472                     IF ( znd < 1.0 ) THEN
1473                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
1474                     ELSE
1475                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
1476                     ENDIF
1477                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
1478                  END DO
1479               ENDIF
1480               !
1481            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1482         END DO
1483      END DO
1484    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
1485
1486    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
1487     !!---------------------------------------------------------------------
1488     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
1489     !!
1490     !! ** Purpose : Calculates the rate at which hbl changes.
1491     !!
1492     !! ** Method  :
1493     !!
1494     !!----------------------------------------------------------------------
1495
1496    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl
1497
1498    INTEGER :: jj, ji
1499    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert
1500    REAL(wp) :: zvel_max, zwb_min
1501    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
1502    REAL(wp) :: zzeta_m = 0.3
1503    REAL(wp) :: zgamma_c = 2.0
1504    REAL(wp) :: zdhoh = 0.1
1505    REAL(wp) :: alpha_bc = 0.5
1506
1507    DO jj = 2, jpjm1
1508       DO ji = 2, jpim1
1509          IF ( lconv(ji,jj) ) THEN    ! Convective
1510             ! Alan is this right?  Yes, it's a bit different from the previous relationship
1511             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) &
1512             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1513             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
1514             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
1515             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
1516             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
1517             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
1518                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
1519
1520             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
1521                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
1522                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
1523                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1524             !
1525             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj)
1526
1527             IF ( ln_osm_mle ) THEN
1528                  !  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 + &
1529                ! &            ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 )          ! Fox-Kemper buoyancy flux average over OSBL
1530                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
1531                   zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  &
1532                        (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
1533                ELSE
1534                   zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
1535                ENDIF
1536                zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1537                IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
1538                   ! OSBL is deepening, entrainment > restratification
1539                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1540                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1541                   ELSE
1542                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15)
1543                   ENDIF
1544                ELSE
1545                   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
1546                   zdhdt(ji,jj) = - zvel_mle(ji,jj)
1547
1548
1549                ENDIF
1550
1551             ELSE
1552                ! Fox-Kemper not used.
1553
1554                  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) / &
1555                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1556                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1557                ! added ajgn 23 July as temporay fix
1558
1559             ENDIF
1560
1561             zdhdt_2(ji,jj) = 0._wp
1562
1563                ! commented out ajgn 23 July as temporay fix
1564!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1565! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary.
1566!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1567!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj)
1568!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1569!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj)
1570!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min )
1571!                     ! Alan no idea what this should be?
1572!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) &
1573!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) &
1574!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj)
1575!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
1576!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN
1577!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj)
1578!                 ENDIF
1579            ELSE                        ! Stable
1580                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
1581                zdhdt_2(ji,jj) = 0._wp
1582                IF ( zdhdt(ji,jj) < 0._wp ) THEN
1583                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
1584                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
1585                ELSE
1586                    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) )
1587                ENDIF
1588                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert
1589            ENDIF
1590         END DO
1591      END DO
1592    END SUBROUTINE zdf_osm_calculate_dhdt
1593
1594    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
1595     !!---------------------------------------------------------------------
1596     !!                   ***  ROUTINE zdf_osm_timestep_hbl  ***
1597     !!
1598     !! ** Purpose : Increments hbl.
1599     !!
1600     !! ** Method  : If thechange in hbl exceeds one model level the change is
1601     !!              is calculated by moving down the grid, changing the buoyancy
1602     !!              jump. This is to ensure that the change in hbl does not
1603     !!              overshoot a stable layer.
1604     !!
1605     !!----------------------------------------------------------------------
1606
1607
1608    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl.
1609
1610    INTEGER :: jk, jj, ji, jm
1611    REAL(wp) :: zhbl_s, zvel_max, zdb
1612    REAL(wp) :: zthermal, zbeta
1613
1614     DO jj = 2, jpjm1
1615         DO ji = 2, jpim1
1616           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
1617!
1618! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
1619!
1620              zhbl_s = hbl(ji,jj)
1621              jm = imld(ji,jj)
1622              zthermal = rab_n(ji,jj,1,jp_tem)
1623              zbeta = rab_n(ji,jj,1,jp_sal)
1624
1625
1626              IF ( lconv(ji,jj) ) THEN
1627!unstable
1628
1629                 IF( ln_osm_mle ) THEN
1630                    zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
1631                 ELSE
1632
1633                    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) / &
1634                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1635
1636                 ENDIF
1637
1638                 DO jk = imld(ji,jj), ibld(ji,jj)
1639                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
1640                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
1641                         &  0.0 ) + zvel_max
1642
1643
1644                    IF ( ln_osm_mle ) THEN
1645                       zhbl_s = zhbl_s + MIN( &
1646                        & 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) ), &
1647                        & e3w_n(ji,jj,jm) )
1648                    ELSE
1649                      zhbl_s = zhbl_s + MIN( &
1650                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
1651                        & e3w_n(ji,jj,jm) )
1652                    ENDIF
1653
1654                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1655
1656                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
1657                 END DO
1658                 hbl(ji,jj) = zhbl_s
1659                 ibld(ji,jj) = jm
1660             ELSE
1661! stable
1662                 DO jk = imld(ji,jj), ibld(ji,jj)
1663                    zdb = MAX( &
1664                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&
1665                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
1666                         & 0.0 ) + &
1667             &       2.0 * zvstr(ji,jj)**2 / zhbl_s
1668
1669                    ! Alan is thuis right? I have simply changed hbli to hbl
1670                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
1671                    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) ) ) * &
1672               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
1673                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
1674                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
1675
1676                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1677                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
1678                 END DO
1679             ENDIF   ! IF ( lconv )
1680             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
1681             ibld(ji,jj) = MAX(jm, 4 )
1682           ELSE
1683! change zero or one model level.
1684             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
1685           ENDIF
1686           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
1687         END DO
1688      END DO
1689
1690    END SUBROUTINE zdf_osm_timestep_hbl
1691
1692    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
1693      !!---------------------------------------------------------------------
1694      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  ***
1695      !!
1696      !! ** Purpose : Calculates thickness of the pycnocline
1697      !!
1698      !! ** Method  : The thickness is calculated from a prognostic equation
1699      !!              that relaxes the pycnocine thickness to a diagnostic
1700      !!              value. The time change is calculated assuming the
1701      !!              thickness relaxes exponentially. This is done to deal
1702      !!              with large timesteps.
1703      !!
1704      !!----------------------------------------------------------------------
1705
1706      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness.
1707      !
1708      INTEGER :: jj, ji
1709      INTEGER :: inhml
1710      REAL(wp) :: zari, ztau, zddhdt
1711
1712
1713      DO jj = 2, jpjm1
1714         DO ji = 2, jpim1
1715            IF ( lconv(ji,jj) ) THEN
1716
1717               IF( ln_osm_mle ) THEN
1718                  IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN
1719                     ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
1720                     IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1721                        IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
1722                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1723                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1724                        ELSE                                                     ! unstable
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 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1727                        ENDIF
1728                        ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
1729                        zddhdt = zari * hbl(ji,jj)
1730                     ELSE
1731                        ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1732                        zddhdt = 0.2 * hbl(ji,jj)
1733                     ENDIF
1734                  ELSE
1735                     ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1736                     zddhdt = 0.2 * hbl(ji,jj)
1737                  ENDIF
1738               ELSE ! ln_osm_mle
1739                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1740                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
1741                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1742                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1743                     ELSE                                                     ! unstable
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 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1746                     ENDIF
1747                     ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
1748                     zddhdt = zari * hbl(ji,jj)
1749                  ELSE
1750                     ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1751                     zddhdt = 0.2 * hbl(ji,jj)
1752                  ENDIF
1753
1754               END IF  ! ln_osm_mle
1755
1756               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1757               ! Alan: this hml is never defined or used
1758            ELSE   ! IF (lconv)
1759               ztau = hbl(ji,jj) / zvstr(ji,jj)
1760               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
1761                  ! boundary layer deepening
1762                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1763                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
1764                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
1765                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
1766                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj)
1767                  ELSE
1768                     zddhdt = 0.2 * hbl(ji,jj)
1769                  ENDIF
1770               ELSE     ! IF(dhdt < 0)
1771                  zddhdt = 0.2 * hbl(ji,jj)
1772               ENDIF    ! IF (dhdt >= 0)
1773               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1774               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
1775               ! Alan: this hml is never defined or used -- do we need it?
1776            ENDIF       ! IF (lconv)
1777
1778            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
1779            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
1780            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
1781            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
1782            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1783         END DO
1784      END DO
1785
1786    END SUBROUTINE zdf_osm_pycnocline_thickness
1787
1788
1789   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )
1790      !!----------------------------------------------------------------------
1791      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  ***
1792      !!
1793      !! ** Purpose :   Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.
1794      !!
1795      !! ** Method  :
1796      !!
1797      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
1798      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
1799
1800
1801      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points
1802      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == !
1803      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy
1804
1805      INTEGER  ::   ji, jj, jk          ! dummy loop indices
1806      INTEGER  ::   ii, ij, ik, ikmax   ! local integers
1807      REAL(wp)                         :: zc
1808      REAL(wp)                         :: zN2_c           ! local buoyancy difference from 10m value
1809      REAL(wp), DIMENSION(jpi,jpj)     :: ztm, zsm, zLf_NH, zLf_MH
1810      REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv
1811      REAL(wp), DIMENSION(jpi,jpj)     :: zmld_midu, zmld_midv
1812!!----------------------------------------------------------------------
1813      !
1814      !                                      !==  MLD used for MLE  ==!
1815
1816      mld_prof(:,:)  = nlb10               ! Initialization to the number of w ocean point
1817      zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
1818      zN2_c = grav * rn_osm_mle_rho_c * r1_rau0   ! convert density criteria into N^2 criteria
1819      DO jk = nlb10, jpkm1
1820         DO jj = 1, jpj                ! Mixed layer level: w-level
1821            DO ji = 1, jpi
1822               ikt = mbkt(ji,jj)
1823               zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
1824               IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
1825            END DO
1826         END DO
1827      END DO
1828      DO jj = 1, jpj
1829         DO ji = 1, jpi
1830            mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
1831            zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
1832         END DO
1833      END DO
1834      ! ensure mld_prof .ge. ibld
1835      !
1836      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation
1837      !
1838      ztm(:,:) = 0._wp
1839      zsm(:,:) = 0._wp
1840      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
1841         DO jj = 1, jpj
1842            DO ji = 1, jpi
1843               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
1844               ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem)
1845               zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal)
1846            END DO
1847         END DO
1848      END DO
1849      ! average temperature and salinity.
1850      ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
1851      zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
1852      ! calculate horizontal gradients at u & v points
1853
1854      DO jj = 2, jpjm1
1855         DO ji = 1, jpim1
1856            zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
1857            zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
1858            zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))
1859            ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )
1860            ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )
1861         END DO
1862      END DO
1863
1864      DO jj = 1, jpjm1
1865         DO ji = 2, jpim1
1866            zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
1867            zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
1868            zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))
1869            ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )
1870            ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )
1871         END DO
1872      END DO
1873
1874      CALL eos_rab(ztsm_midu, zmld_midu, zabu)
1875      CALL eos_rab(ztsm_midv, zmld_midv, zabv)
1876
1877      DO jj = 2, jpjm1
1878         DO ji = 1, jpim1
1879            dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
1880         END DO
1881      END DO
1882      DO jj = 1, jpjm1
1883         DO ji = 2, jpim1
1884            dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
1885         END DO
1886      END DO
1887
1888 END SUBROUTINE zdf_osm_zmld_horizontal_gradients
1889  SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )
1890      !!----------------------------------------------------------------------
1891      !!                  ***  ROUTINE zdf_osm_mle_parameters  ***
1892      !!
1893      !! ** Purpose :   Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.
1894      !!
1895      !! ** Method  :
1896      !!
1897      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
1898      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
1899
1900      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zwb_fk, zvel_mle, zdiff_mle
1901      INTEGER  ::   ji, jj, jk          ! dummy loop indices
1902      INTEGER  ::   ii, ij, ik  ! local integers
1903      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
1904      REAL(wp) ::   zdb_mle, ztmp, zdbds_mle
1905
1906      mld_prof(:,:) = 4
1907      DO jk = 5, jpkm1
1908        DO jj = 2, jpjm1
1909          DO ji = 2, jpim1
1910            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
1911          END DO
1912        END DO
1913      END DO
1914      ! DO jj = 2, jpjm1
1915      !    DO ji = 1, jpim1
1916      !       zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
1917      !    END DO
1918      !  END DO
1919   ! Timestep mixed layer eddy depth.
1920      DO jj = 2, jpjm1
1921        DO ji = 2, jpim1
1922          zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG
1923          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
1924             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.
1925          ELSE
1926            ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10
1927             ! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN
1928             !   hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau
1929             ! ELSE
1930             !   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
1931             ! ENDIF
1932             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
1933               hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau
1934             ELSE
1935               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
1936             ENDIF
1937          ENDIF
1938          hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj))
1939        END DO
1940      END DO
1941
1942      mld_prof = 4
1943      DO jk = 5, jpkm1
1944        DO jj = 2, jpjm1
1945          DO ji = 2, jpim1
1946            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
1947          END DO
1948        END DO
1949      END DO
1950      DO jj = 2, jpjm1
1951         DO ji = 2, jpim1
1952            zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj))
1953         END DO
1954       END DO
1955   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
1956
1957      DO jj = 2, jpjm1
1958        DO ji = 2, jpim1
1959          IF ( lconv(ji,jj) ) THEN
1960             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
1961             ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) &
1962             !      & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) )
1963             ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) &
1964             !      & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) )
1965             zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &
1966                  & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
1967             zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle
1968      ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
1969             zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
1970             zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf
1971          ENDIF
1972        END DO
1973      END DO
1974END SUBROUTINE zdf_osm_mle_parameters
1975
1976END SUBROUTINE zdf_osm
1977
1978
1979   SUBROUTINE zdf_osm_init
1980     !!----------------------------------------------------------------------
1981     !!                  ***  ROUTINE zdf_osm_init  ***
1982     !!
1983     !! ** Purpose :   Initialization of the vertical eddy diffivity and
1984     !!      viscosity when using a osm turbulent closure scheme
1985     !!
1986     !! ** Method  :   Read the namosm namelist and check the parameters
1987     !!      called at the first timestep (nit000)
1988     !!
1989     !! ** input   :   Namlist namosm
1990     !!----------------------------------------------------------------------
1991     INTEGER  ::   ios            ! local integer
1992     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1993     REAL z1_t2
1994     !!
1995     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
1996          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
1997          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
1998          & ,ln_osm_mle
1999! Namelist for Fox-Kemper parametrization.
2000      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
2001           & rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau
2002
2003     !!----------------------------------------------------------------------
2004     !
2005     REWIND( numnam_ref )              ! Namelist namzdf_osm in reference namelist : Osmosis ML model
2006     READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
2007901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist', lwp )
2008
2009     REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
2010     READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
2011902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist', lwp )
2012     IF(lwm) WRITE ( numond, namzdf_osm )
2013
2014     IF(lwp) THEN                    ! Control print
2015        WRITE(numout,*)
2016        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
2017        WRITE(numout,*) '~~~~~~~~~~~~'
2018        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters'
2019        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la
2020        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle
2021        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
2022        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
2023        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes
2024        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
2025        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave
2026        SELECT CASE (nn_osm_wave)
2027        CASE(0)
2028           WRITE(numout,*) '     calculated assuming constant La#=0.3'
2029        CASE(1)
2030           WRITE(numout,*) '     calculated from Pierson Moskowitz wind-waves'
2031        CASE(2)
2032           WRITE(numout,*) '     calculated from ECMWF wave fields'
2033        END SELECT
2034        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
2035        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix
2036        WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
2037        WRITE(numout,*) '     maximum shear diffusivity at Rig = 0    (m2/s) rn_difri = ', rn_difri
2038        WRITE(numout,*) '     Use large mixing below BL when unstable       ln_convmix = ', ln_convmix
2039        WRITE(numout,*) '     diffusivity when unstable below BL     (m2/s) rn_difconv = ', rn_difconv
2040     ENDIF
2041
2042     !                              ! allocate zdfosm arrays
2043     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
2044
2045
2046     IF( ln_osm_mle ) THEN
2047! Initialise Fox-Kemper parametrization
2048         REWIND( numnam_ref )              ! Namelist namosm_mle in reference namelist : Tracer advection scheme
2049         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
2050903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist',lwp)
2051
2052         REWIND( numnam_cfg )              ! Namelist namosm_mle in configuration namelist : Tracer advection scheme
2053         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
2054904      IF( ios >  0 )   CALL ctl_nam ( ios , 'namosm_mle in configuration namelist',lwp)
2055         IF(lwm) WRITE ( numond, namosm_mle )
2056
2057         IF(lwp) THEN                     ! Namelist print
2058            WRITE(numout,*)
2059            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
2060            WRITE(numout,*) '~~~~~~~~~~~~~'
2061            WRITE(numout,*) '   Namelist namosm_mle : '
2062            WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_osm_mle    = ', nn_osm_mle
2063            WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_osm_mle_ce    = ', rn_osm_mle_ce
2064            WRITE(numout,*) '         scale of ML front (ML radius of deformation) (nn_osm_mle=0)      rn_osm_mle_lf     = ', rn_osm_mle_lf, 'm'
2065            WRITE(numout,*) '         maximum time scale of MLE                    (nn_osm_mle=0)      rn_osm_mle_time   = ', rn_osm_mle_time, 's'
2066            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg'
2067            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c
2068            WRITE(numout,*) '         Threshold used to define ML for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s'
2069            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's'
2070         ENDIF         !
2071     ENDIF
2072      !
2073      IF(lwp) THEN
2074         WRITE(numout,*)
2075         IF( ln_osm_mle ) THEN
2076            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
2077            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
2078            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation'
2079         ELSE
2080            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
2081         ENDIF
2082      ENDIF
2083      !
2084      IF( ln_osm_mle ) THEN                ! MLE initialisation
2085         !
2086         rb_c = grav * rn_osm_mle_rho_c /rau0        ! Mixed Layer buoyancy criteria
2087         IF(lwp) WRITE(numout,*)
2088         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
2089         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
2090         !
2091         IF( nn_osm_mle == 0 ) THEN           ! MLE array allocation & initialisation            !
2092!
2093         ELSEIF( nn_osm_mle == 1 ) THEN           ! MLE array allocation & initialisation
2094            rc_f = rn_osm_mle_ce/ (  5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat )  )
2095            !
2096         ENDIF
2097         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
2098         z1_t2 = 2.e-5
2099         do jj=1,jpj
2100            do ji = 1,jpi
2101               r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
2102            end do
2103         end do
2104         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
2105         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
2106         !
2107      ENDIF
2108
2109     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle
2110
2111
2112     IF( ln_zdfddm) THEN
2113        IF(lwp) THEN
2114           WRITE(numout,*)
2115           WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
2116           WRITE(numout,*) '    CAUTION : done in routine zdfosm, not in routine zdfddm '
2117        ENDIF
2118     ENDIF
2119
2120
2121     !set constants not in namelist
2122     !-----------------------------
2123
2124     IF(lwp) THEN
2125        WRITE(numout,*)
2126     ENDIF
2127
2128     IF (nn_osm_wave == 0) THEN
2129        dstokes(:,:) = rn_osm_dstokes
2130     END IF
2131
2132     ! Horizontal average : initialization of weighting arrays
2133     ! -------------------
2134
2135     SELECT CASE ( nn_ave )
2136
2137     CASE ( 0 )                ! no horizontal average
2138        IF(lwp) WRITE(numout,*) '          no horizontal average on avt'
2139        IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
2140        ! weighting mean arrays etmean
2141        !           ( 1  1 )
2142        ! avt = 1/4 ( 1  1 )
2143        !
2144        etmean(:,:,:) = 0.e0
2145
2146        DO jk = 1, jpkm1
2147           DO jj = 2, jpjm1
2148              DO ji = 2, jpim1   ! vector opt.
2149                 etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
2150                      &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
2151                      &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
2152              END DO
2153           END DO
2154        END DO
2155
2156     CASE ( 1 )                ! horizontal average
2157        IF(lwp) WRITE(numout,*) '          horizontal average on avt'
2158        ! weighting mean arrays etmean
2159        !           ( 1/2  1  1/2 )
2160        ! avt = 1/8 ( 1    2  1   )
2161        !           ( 1/2  1  1/2 )
2162        etmean(:,:,:) = 0.e0
2163
2164        DO jk = 1, jpkm1
2165           DO jj = 2, jpjm1
2166              DO ji = 2, jpim1   ! vector opt.
2167                 etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
2168                      & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
2169                      &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
2170                      &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
2171                      &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
2172                      &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
2173              END DO
2174           END DO
2175        END DO
2176
2177     CASE DEFAULT
2178        WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
2179        CALL ctl_stop( ctmp1 )
2180
2181     END SELECT
2182
2183     ! Initialization of vertical eddy coef. to the background value
2184     ! -------------------------------------------------------------
2185     DO jk = 1, jpk
2186        avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
2187     END DO
2188
2189     ! zero the surface flux for non local term and osm mixed layer depth
2190     ! ------------------------------------------------------------------
2191     ghamt(:,:,:) = 0.
2192     ghams(:,:,:) = 0.
2193     ghamu(:,:,:) = 0.
2194     ghamv(:,:,:) = 0.
2195     !
2196     IF( lwxios ) THEN
2197        CALL iom_set_rstw_var_active('wn')
2198        CALL iom_set_rstw_var_active('hbl')
2199        CALL iom_set_rstw_var_active('dh')
2200        IF( ln_osm_mle ) THEN
2201            CALL iom_set_rstw_var_active('hmle')
2202        END IF
2203     ENDIF
2204   END SUBROUTINE zdf_osm_init
2205
2206
2207   SUBROUTINE osm_rst( kt, cdrw )
2208     !!---------------------------------------------------------------------
2209     !!                   ***  ROUTINE osm_rst  ***
2210     !!
2211     !! ** Purpose :   Read or write BL fields in restart file
2212     !!
2213     !! ** Method  :   use of IOM library. If the restart does not contain
2214     !!                required fields, they are recomputed from stratification
2215     !!----------------------------------------------------------------------
2216
2217     INTEGER, INTENT(in) :: kt
2218     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
2219
2220     INTEGER ::   id1, id2, id3   ! iom enquiry index
2221     INTEGER  ::   ji, jj, jk     ! dummy loop indices
2222     INTEGER  ::   iiki, ikt ! local integer
2223     REAL(wp) ::   zhbf           ! tempory scalars
2224     REAL(wp) ::   zN2_c           ! local scalar
2225     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
2226     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
2227     !!----------------------------------------------------------------------
2228     !
2229     !!-----------------------------------------------------------------------------
2230     ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
2231     !!-----------------------------------------------------------------------------
2232     IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
2233        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. )
2234        IF( id1 > 0 ) THEN                       ! 'wn' exists; read
2235           CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios )
2236           WRITE(numout,*) ' ===>>>> :  wn read from restart file'
2237        ELSE
2238           wn(:,:,:) = 0._wp
2239           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
2240        END IF
2241
2242        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
2243        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. )
2244        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return
2245           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
2246           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  )
2247           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file'
2248           IF( ln_osm_mle ) THEN
2249              id3 = iom_varid( numror, 'hmle'   , ldstop = .FALSE. )
2250              IF( id3 > 0) THEN
2251                 CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios )
2252                 WRITE(numout,*) ' ===>>>> :  hmle read from restart file'
2253              ELSE
2254                 WRITE(numout,*) ' ===>>>> :  hmle not found, set to hbl'
2255                 hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
2256              END IF
2257           END IF
2258           RETURN
2259        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate
2260           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
2261        END IF
2262     END IF
2263
2264     !!-----------------------------------------------------------------------------
2265     ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
2266     !!-----------------------------------------------------------------------------
2267     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
2268        IF(lwp) WRITE(numout,*) '---- osm-rst ----'
2269         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn,   ldxios = lwxios )
2270         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl,  ldxios = lwxios )
2271         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh,   ldxios = lwxios )
2272         IF( ln_osm_mle ) THEN
2273            CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios )
2274         END IF
2275        RETURN
2276     END IF
2277
2278     !!-----------------------------------------------------------------------------
2279     ! Getting hbl, no restart file with hbl, so calculate from surface stratification
2280     !!-----------------------------------------------------------------------------
2281     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
2282     ! w-level of the mixing and mixed layers
2283     CALL eos_rab( tsn, rab_n )
2284     CALL bn2(tsn, rab_n, rn2)
2285     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
2286     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
2287     zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
2288     !
2289     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
2290     DO jk = 1, jpkm1
2291        DO jj = 1, jpj              ! Mixed layer level: w-level
2292           DO ji = 1, jpi
2293              ikt = mbkt(ji,jj)
2294              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
2295              IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
2296           END DO
2297        END DO
2298     END DO
2299     !
2300     DO jj = 1, jpj
2301        DO ji = 1, jpi
2302           iiki = MAX(4,imld_rst(ji,jj))
2303           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth
2304           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth
2305        END DO
2306     END DO
2307
2308     WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
2309
2310     IF( ln_osm_mle ) THEN
2311        hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
2312        WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
2313     END IF
2314
2315     wn(:,:,:) = 0._wp
2316     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
2317   END SUBROUTINE osm_rst
2318
2319
2320   SUBROUTINE tra_osm( kt )
2321      !!----------------------------------------------------------------------
2322      !!                  ***  ROUTINE tra_osm  ***
2323      !!
2324      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
2325      !!
2326      !! ** Method  :   ???
2327      !!----------------------------------------------------------------------
2328      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
2329      !!----------------------------------------------------------------------
2330      INTEGER, INTENT(in) :: kt
2331      INTEGER :: ji, jj, jk
2332      !
2333      IF( kt == nit000 ) THEN
2334         IF(lwp) WRITE(numout,*)
2335         IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
2336         IF(lwp) WRITE(numout,*) '~~~~~~~   '
2337      ENDIF
2338
2339      IF( l_trdtra )   THEN                    !* Save ta and sa trends
2340         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
2341         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
2342      ENDIF
2343
2344      DO jk = 1, jpkm1
2345         DO jj = 2, jpjm1
2346            DO ji = 2, jpim1
2347               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
2348                  &                 - (  ghamt(ji,jj,jk  )  &
2349                  &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
2350               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
2351                  &                 - (  ghams(ji,jj,jk  )  &
2352                  &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
2353            END DO
2354         END DO
2355      END DO
2356
2357      ! save the non-local tracer flux trends for diagnostics
2358      IF( l_trdtra )   THEN
2359         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
2360         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
2361
2362         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt )
2363         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds )
2364         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
2365      ENDIF
2366
2367      IF(ln_ctl) THEN
2368         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   &
2369         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
2370      ENDIF
2371      !
2372   END SUBROUTINE tra_osm
2373
2374
2375   SUBROUTINE trc_osm( kt )          ! Dummy routine
2376      !!----------------------------------------------------------------------
2377      !!                  ***  ROUTINE trc_osm  ***
2378      !!
2379      !! ** Purpose :   compute and add to the passive tracer trend the non-local
2380      !!                 passive tracer flux
2381      !!
2382      !!
2383      !! ** Method  :   ???
2384      !!----------------------------------------------------------------------
2385      !
2386      !!----------------------------------------------------------------------
2387      INTEGER, INTENT(in) :: kt
2388      WRITE(*,*) 'trc_osm: Not written yet', kt
2389   END SUBROUTINE trc_osm
2390
2391
2392   SUBROUTINE dyn_osm( kt )
2393      !!----------------------------------------------------------------------
2394      !!                  ***  ROUTINE dyn_osm  ***
2395      !!
2396      !! ** Purpose :   compute and add to the velocity trend the non-local flux
2397      !! copied/modified from tra_osm
2398      !!
2399      !! ** Method  :   ???
2400      !!----------------------------------------------------------------------
2401      INTEGER, INTENT(in) ::   kt   !
2402      !
2403      INTEGER :: ji, jj, jk   ! dummy loop indices
2404      !!----------------------------------------------------------------------
2405      !
2406      IF( kt == nit000 ) THEN
2407         IF(lwp) WRITE(numout,*)
2408         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
2409         IF(lwp) WRITE(numout,*) '~~~~~~~   '
2410      ENDIF
2411      !code saving tracer trends removed, replace with trdmxl_oce
2412
2413      DO jk = 1, jpkm1           ! add non-local u and v fluxes
2414         DO jj = 2, jpjm1
2415            DO ji = 2, jpim1
2416               ua(ji,jj,jk) =  ua(ji,jj,jk)                      &
2417                  &                 - (  ghamu(ji,jj,jk  )  &
2418                  &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
2419               va(ji,jj,jk) =  va(ji,jj,jk)                      &
2420                  &                 - (  ghamv(ji,jj,jk  )  &
2421                  &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
2422            END DO
2423         END DO
2424      END DO
2425      !
2426      ! code for saving tracer trends removed
2427      !
2428   END SUBROUTINE dyn_osm
2429
2430   !!======================================================================
2431
2432END MODULE zdfosm
Note: See TracBrowser for help on using the repository browser.