New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfosm.F90 in NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF – NEMO

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

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

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