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 @ 14409

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

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