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/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90 @ 14408

Last change on this file since 14408 was 14408, checked in by agn, 3 years ago

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