source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90

Last change on this file was 12928, checked in by smueller, 5 months ago

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

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