source: NEMO/branches/UKMO/NEMO_4.0_OSMOSIS/src/OCE/ZDF/zdfosm.F90 @ 11860

Last change on this file since 11860 was 11860, checked in by cguiavarch, 16 months ago

Update with an updated version of zdfosm.F90 from George Nurser

File size: 103.8 KB
Line 
1MODULE zdfosm
2   !!======================================================================
3   !!                       ***  MODULE  zdfosm  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the OSMOSIS
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !!  History : NEMO 4.0  ! A. Grant, G. Nurser
8   !! 15/03/2017  Changed calculation of pycnocline thickness in unstable conditions and stable conditions AG
9   !! 15/03/2017  Calculation of pycnocline gradients for stable conditions changed. Pycnocline gradients now depend on stability of the OSBL. A.G
10   !! 06/06/2017  (1) Checks on sign of buoyancy jump in calculation of  OSBL depth.  A.G.
11   !!             (2) Removed variable zbrad0, zbradh and zbradav since they are not used.
12   !!             (3) Approximate treatment for shear turbulence.
13   !!                        Minimum values for zustar and zustke.
14   !!                        Add velocity scale, zvstr, that tends to zustar for large Langmuir numbers.
15   !!                        Limit maximum value for Langmuir number.
16   !!                        Use zvstr in definition of stability parameter zhol.
17   !!             (4) Modified parametrization of entrainment flux, changing original coefficient 0.0485 for Langmuir contribution to 0.135 * zla
18   !!             (5) For stable boundary layer add factor that depends on length of timestep to 'slow' collapse and growth. Make sure buoyancy jump not negative.
19   !!             (6) For unstable conditions when growth is over multiple levels, limit change to maximum of one level per cycle through loop.
20   !!             (7) Change lower limits for loops that calculate OSBL averages from 1 to 2. Large gradients between levels 1 and 2 can cause problems.
21   !!             (8) Change upper limits from ibld-1 to ibld.
22   !!             (9) Calculation of pycnocline thickness in unstable conditions. Check added to ensure that buoyancy jump is positive before calculating Ri.
23   !!            (10) Thickness of interface layer at base of the stable OSBL set by Richardson number. Gives continuity in transition from unstable OSBL.
24   !!            (11) Checks that buoyancy jump is poitive when calculating pycnocline profiles.
25   !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity.
26   !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information
27   !!            (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence.
28   !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added.
29   !!            (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out)
30   !!            (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out)
31   !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made,
32   !!                  (a) Pycnocline temperature and salinity profies changed for unstable layers
33   !!                  (b) The stable OSBL depth parametrization changed.
34   !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code.
35   !! 23/05/19   (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1
36   !!----------------------------------------------------------------------
37
38   !!----------------------------------------------------------------------
39   !!   'ln_zdfosm'                                             OSMOSIS scheme
40   !!----------------------------------------------------------------------
41   !!   zdf_osm       : update momentum and tracer Kz from osm scheme
42   !!   zdf_osm_init  : initialization, namelist read, and parameters control
43   !!   osm_rst       : read (or initialize) and write osmosis restart fields
44   !!   tra_osm       : compute and add to the T & S trend the non-local flux
45   !!   trc_osm       : compute and add to the passive tracer trend the non-local flux (TBD)
46   !!   dyn_osm       : compute and add to u & v trensd the non-local flux
47   !!
48   !! Subroutines in revised code.
49   !!----------------------------------------------------------------------
50   USE oce            ! ocean dynamics and active tracers
51                      ! uses 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   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux
82   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv    !: non-local v-momentum flux
83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamt    !: non-local temperature flux (gamma/<ws>o)
84   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghams    !: non-local salinity flux (gamma/<ws>o)
85   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt
86   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth
87   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh   ! depth of pycnocline
88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml  ! ML depth
89   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift.
90   !                      !!** Namelist  namzdf_osm  **
91   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la
92   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number
93   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift
94   REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs
95   INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt
96   INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
97   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la
98
99
100   LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing
101   REAL(wp) ::   rn_riinfty   = 0.7     ! local Richardson Number limit for shear instability
102   REAL(wp) ::   rn_difri    =  0.005   ! maximum shear mixing at Rig = 0    (m2/s)
103   LOGICAL  ::   ln_convmix  = .true.   ! Convective instability mixing
104   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s)
105
106   !                                    !!! ** General constants  **
107   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero
108   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw
109   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3
110   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3
111
112   INTEGER :: idebug = 236
113   INTEGER :: jdebug = 228
114   !!----------------------------------------------------------------------
115   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
116   !! $Id$
117   !! Software governed by the CeCILL license (see ./LICENSE)
118   !!----------------------------------------------------------------------
119CONTAINS
120
121   INTEGER FUNCTION zdf_osm_alloc()
122      !!----------------------------------------------------------------------
123      !!                 ***  FUNCTION zdf_osm_alloc  ***
124      !!----------------------------------------------------------------------
125     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &
126          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
127          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
128!     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &    ! would ths be better ?
129!          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
130!          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
131!     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
132!     IF ( ln_osm_mle ) THEN
133!        Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc )
134!        IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays')
135!     ENDIF
136
137     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
138     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
139   END FUNCTION zdf_osm_alloc
140
141
142   SUBROUTINE zdf_osm( kt, p_avm, p_avt )
143      !!----------------------------------------------------------------------
144      !!                   ***  ROUTINE zdf_osm  ***
145      !!
146      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
147      !!      coefficients and non local mixing using the OSMOSIS scheme
148      !!
149      !! ** Method :   The boundary layer depth hosm is diagnosed at tracer points
150      !!      from profiles of buoyancy, and shear, and the surface forcing.
151      !!      Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
152      !!
153      !!                      Kx =  hosm  Wx(sigma) G(sigma)
154      !!
155      !!             and the non local term ghamt = Cs / Ws(sigma) / hosm
156      !!      Below hosm  the coefficients are the sum of mixing due to internal waves
157      !!      shear instability and double diffusion.
158      !!
159      !!      -1- Compute the now interior vertical mixing coefficients at all depths.
160      !!      -2- Diagnose the boundary layer depth.
161      !!      -3- Compute the now boundary layer vertical mixing coefficients.
162      !!      -4- Compute the now vertical eddy vicosity and diffusivity.
163      !!      -5- Smoothing
164      !!
165      !!        N.B. The computation is done from jk=2 to jpkm1
166      !!             Surface value of avt are set once a time to zero
167      !!             in routine zdf_osm_init.
168      !!
169      !! ** Action  :   update the non-local terms ghamts
170      !!                update avt (before vertical eddy coef.)
171      !!
172      !! References : Large W.G., Mc Williams J.C. and Doney S.C.
173      !!         Reviews of Geophysics, 32, 4, November 1994
174      !!         Comments in the code refer to this paper, particularly
175      !!         the equation number. (LMD94, here after)
176      !!----------------------------------------------------------------------
177      INTEGER                   , INTENT(in   ) ::   kt            ! ocean time step
178      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points)
179      !!
180      INTEGER ::   ji, jj, jk                   ! dummy loop indices
181      INTEGER ::   ikbot, jkmax, jkm1, jkp2     !
182
183      REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube      !
184      REAL(wp) ::   zbeta, zthermal                                  !
185      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
186      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      !
187
188      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density
189      INTEGER  ::   jm                          ! dummy loop indices
190      REAL(wp) ::   zr1, zr2, zr3, zr4, zrhop   ! Compression terms
191      REAL(wp) ::   zflag, zrn2, zdep21, zdep32, zdep43
192      REAL(wp) ::   zesh2, zri, zfri            ! Interior richardson mixing
193      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
194      REAL(wp) :: zt,zs,zu,zv,zrh               ! variables used in constructing averages
195! Scales
196      REAL(wp), DIMENSION(jpi,jpj) :: zrad0     ! Surface solar temperature flux (deg m/s)
197      REAL(wp), DIMENSION(jpi,jpj) :: zradh     ! Radiative flux at bl base (Buoyancy units)
198      REAL(wp), DIMENSION(jpi,jpj) :: zradav    ! Radiative flux, bl average (Buoyancy Units)
199      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
200      REAL(wp), DIMENSION(jpi,jpj) :: zwstrl    ! Langmuir velocity scale
201      REAL(wp), DIMENSION(jpi,jpj) :: zvstr     ! Velocity scale that ends to zustar for large Langmuir number.
202      REAL(wp), DIMENSION(jpi,jpj) :: zwstrc    ! Convective velocity scale
203      REAL(wp), DIMENSION(jpi,jpj) :: zuw0      ! Surface u-momentum flux
204      REAL(wp), DIMENSION(jpi,jpj) :: zvw0      ! Surface v-momentum flux
205      REAL(wp), DIMENSION(jpi,jpj) :: zwth0     ! Surface heat flux (Kinematic)
206      REAL(wp), DIMENSION(jpi,jpj) :: zws0      ! Surface freshwater flux
207      REAL(wp), DIMENSION(jpi,jpj) :: zwb0      ! Surface buoyancy flux
208      REAL(wp), DIMENSION(jpi,jpj) :: zwthav    ! Heat flux - bl average
209      REAL(wp), DIMENSION(jpi,jpj) :: zwsav     ! freshwater flux - bl average
210      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average
211      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux
212
213      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift
214      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number
215      REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
216      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
217      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer
218      LOGICAL, DIMENSION(jpi,jpj)  :: lconv     ! unstable/stable bl
219
220      ! mixed-layer variables
221
222      INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
223      INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
224      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
225      REAL(wp) :: zugrad,zvgrad        ! temporary variables for calculating pycnocline shear
226
227      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid
228      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid
229      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid
230      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
231      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2                                    ! correction to dhdt due to internal structure.
232      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext              ! external temperature/salinity and buoyancy gradients
233      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer
234      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer
235      REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer
236      REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
237      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
238      REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline
239      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline
240      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline
241      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline
242      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc    ! u-shear across the pycnocline
243      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc    ! v-shear across the pycnocline
244
245      ! Flux-gradient relationship variables
246
247      REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale.
248
249      REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity
250      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
251      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.
252      REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
253
254      ! For calculating Ri#-dependent mixing
255      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du   ! u-shear^2
256      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv   ! v-shear^2
257      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
258
259      ! Temporary variables
260      INTEGER :: inhml
261      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
262      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables
263      REAL(wp) :: zthick, zz0, zz1 ! temporary variables
264      REAL(wp) :: zvel_max, zhbl_s ! temporary variables
265      REAL(wp) :: zfac             ! temporary variable
266      REAL(wp) :: zus_x, zus_y     ! temporary Stokes drift
267      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
268      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
269
270      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme
271      REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt
272      REAL(wp) :: zzeta_s = 0._wp
273      REAL(wp) :: zzeta_v = 0.46
274      REAL(wp) :: zabsstke
275
276      ! For debugging
277      INTEGER :: ikt
278      !!--------------------------------------------------------------------
279      !
280      ibld(:,:)   = 0     ; imld(:,:)  = 0
281      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp
282      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp
283      zvw0(:,:)   = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:)      = 0._wp ; zwb0(:,:)      = 0._wp
284      zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:)     = 0._wp ; zwb_ent(:,:)   = 0._wp
285      zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
286      zhol(:,:)   = 0._wp
287      lconv(:,:)  = .FALSE.
288      ! mixed layer
289      ! no initialization of zhbl or zhml (or zdh?)
290      zhbl(:,:)    = 1._wp ; zhml(:,:)    = 1._wp ; zdh(:,:)      = 1._wp ; zdhdt(:,:)   = 0._wp
291      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp
292      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp
293
294      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp
295      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp
296      zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp
297      zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:)  = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp
298      zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp
299      !
300      zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
301      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
302      !
303      zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp
304      ! Flux-Gradient arrays.
305      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp
306      zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp
307      zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp
308      zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp
309      zhbl_t(:,:)     = 0._wp ; zdhdt(:,:)      = 0._wp
310
311      zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
312      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp
313
314      zdhdt_2(:,:) = 0._wp
315      ! hbl = MAX(hbl,epsln)
316      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
317      ! Calculate boundary layer scales
318      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
319
320      ! Assume two-band radiation model for depth of OSBL
321     zz0 =       rn_abs       ! surface equi-partition in 2-bands
322     zz1 =  1. - rn_abs
323     DO jj = 2, jpjm1
324        DO ji = 2, jpim1
325           ! Surface downward irradiance (so always +ve)
326           zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp
327           ! Downwards irradiance at base of boundary layer
328           zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
329           ! Downwards irradiance averaged over depth of the OSBL
330           zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
331                 &                         + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
332        END DO
333     END DO
334     ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
335     DO jj = 2, jpjm1
336        DO ji = 2, jpim1
337           zthermal = rab_n(ji,jj,1,jp_tem)
338           zbeta    = rab_n(ji,jj,1,jp_sal)
339           ! Upwards surface Temperature flux for non-local term
340           zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
341           ! Upwards surface salinity flux for non-local term
342           zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
343           ! Non radiative upwards surface buoyancy flux
344           zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj)
345           ! turbulent heat flux averaged over depth of OSBL
346           zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
347           ! turbulent salinity flux averaged over depth of the OBSL
348           zwsav(ji,jj) = 0.5 * zws0(ji,jj)
349           ! turbulent buoyancy flux averaged over the depth of the OBSBL
350           zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj)
351           ! Surface upward velocity fluxes
352           zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
353           zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
354           ! Friction velocity (zustar), at T-point : LMD94 eq. 2
355           zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
356           zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
357           zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
358        END DO
359     END DO
360     ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
361     SELECT CASE (nn_osm_wave)
362     ! Assume constant La#=0.3
363     CASE(0)
364        DO jj = 2, jpjm1
365           DO ji = 2, jpim1
366              zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
367              zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
368              zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 )
369              dstokes(ji,jj) = rn_osm_dstokes
370              ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
371           END DO
372        END DO
373     ! Assume Pierson-Moskovitz wind-wave spectrum
374     CASE(1)
375        DO jj = 2, jpjm1
376           DO ji = 2, jpim1
377              ! Use wind speed wndm included in sbc_oce module
378              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
379              dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
380           END DO
381        END DO
382     ! Use ECMWF wave fields as output from SBCWAVE
383     CASE(2)
384        zfac =  2.0_wp * rpi / 16.0_wp
385
386        DO jj = 2, jpjm1
387           DO ji = 2, jpim1
388              ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
389              !    The coefficient 0.8 gives La=0.3  in this situation.
390              ! It could represent the effects of the spread of wave directions
391              ! around the mean wind. The effect of this adjustment needs to be tested.
392              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)
393              zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8)
394              dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
395           END DO
396        END DO
397     END SELECT
398
399     ! Langmuir velocity scale (zwstrl), La # (zla)
400     ! mixed scale (zvstr), convective velocity scale (zwstrc)
401     DO jj = 2, jpjm1
402        DO ji = 2, jpim1
403           ! Langmuir velocity scale (zwstrl), at T-point
404           zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
405           zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)
406           IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))
407           ! Velocity scale that tends to zustar for large Langmuir numbers
408           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + &
409                & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
410
411           ! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
412           ! Note zustke and zwstrl are not amended.
413           !
414           ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
415           IF ( zwbav(ji,jj) > 0.0) THEN
416              zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
417              zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
418              lconv(ji,jj) = .TRUE.
419           ELSE
420              zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln )
421              lconv(ji,jj) = .FALSE.
422           ENDIF
423        END DO
424     END DO
425
426     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
427     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
428     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
429     ! BL must be always 4 levels deep.
430      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) )
431      ibld(:,:) = 4
432      DO jk = 5, jpkm1
433         DO jj = 2, jpjm1
434            DO ji = 2, jpim1
435               IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
436                  ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
437               ENDIF
438            END DO
439         END DO
440      END DO
441
442      DO jj = 2, jpjm1
443         DO ji = 2, jpim1
444            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
445            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 ))
446            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
447         END DO
448      END DO
449      ! Averages over well-mixed and boundary layer
450      ! Alan: do we need zb_nl?, zb_ml?
451      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl)
452      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
453! External gradient
454      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )
455
456! Rate of change of hbl
457      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
458      ! Calculate averages over depth of boundary layer
459      imld = ibld           ! use imld to hold previous blayer index
460      ibld(:,:) = 4
461
462         DO jj = 2, jpjm1
463            DO ji = 2, jpim1
464               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
465               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:))
466               !zdhdt(ji,jj) = MIN(zdhdt(ji,jj), (zhbl_t(ji,jj) - hbl(ji,jj))/rn_rdt + wn(ji,jj,ibld(ji,jj)))
467            END DO
468         END DO
469      ! adjustment to represent limiting by ocean bottom
470
471      DO jk = 4, jpkm1
472         DO jj = 2, jpjm1
473            DO ji = 2, jpim1
474               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
475                  ibld(ji,jj) = jk
476               ENDIF
477            END DO
478         END DO
479      END DO
480
481!
482! Step through model levels taking account of buoyancy change to determine the effect on dhdt
483!
484      CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
485      ! Alan: do we need zb_ml?
486      CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )
487!
488!
489      CALL zdf_osm_pycnocline_thickness( dh, zdh )
490      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms.
491!
492    ! Average over the depth of the mixed layer in the convective boundary layer
493    ! Alan: do we need zb_ml?
494    CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )
495    ! rotate mean currents and changes onto wind align co-ordinates
496    !
497    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
498    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
499     zuw_bse = 0._wp
500     zvw_bse = 0._wp
501     zwth_ent = 0._wp
502     zws_ent = 0._wp
503     DO jj = 2, jpjm1
504        DO ji = 2, jpim1
505           IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN
506              IF ( lconv(ji,jj) ) THEN
507             zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + &
508                      &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ &
509                      &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln)
510            zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ &
511                      &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj))
512                 IF ( zdb_ml(ji,jj) > 0._wp ) THEN
513                    zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
514                    zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
515                 ENDIF
516              ELSE
517                 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )
518                 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )
519              ENDIF
520           ENDIF
521        END DO
522     END DO
523
524      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
525      !  Pycnocline gradients for scalars and velocity
526      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
527
528      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )
529      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc )
530      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )
531       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
532       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
533       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
534
535       DO jj = 2, jpjm1
536          DO ji = 2, jpim1
537             IF ( lconv(ji,jj) ) THEN
538               zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
539               zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
540               zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
541               zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
542               zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third
543               zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln )
544             ELSE
545               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
546               zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
547             END IF
548          END DO
549       END DO
550!
551       DO jj = 2, jpjm1
552          DO ji = 2, jpim1
553             IF ( lconv(ji,jj) ) THEN
554                DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
555                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
556                    !
557                    zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5
558                    !
559                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) &
560                         &            *                                      ( 1.0 -               0.5 * zznd_ml**2 )
561                END DO
562                ! pycnocline - if present linear profile
563                IF ( zdh(ji,jj) > 0._wp ) THEN
564                   zgamma_b = 6.0
565                   DO jk = imld(ji,jj)+1 , ibld(ji,jj)
566                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
567                       !
568                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
569                       !
570                       zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
571                   END DO
572                   IF ( ibld_ext == 0 ) THEN
573                       zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
574                       zviscos(ji,jj,ibld(ji,jj)) = 0._wp
575                   ELSE
576                       zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )
577                       zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )
578                   ENDIF
579                ENDIF
580                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out
581                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj))
582                zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj))
583                ! could be taken out, take account of entrainment represents as a diffusivity
584                ! should remove w from here, represents entrainment
585             ELSE
586             ! stable conditions
587                DO jk = 2, ibld(ji,jj)
588                   zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
589                   zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
590                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
591                END DO
592
593                IF ( ibld_ext == 0 ) THEN
594                   zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
595                   zviscos(ji,jj,ibld(ji,jj)) = 0._wp
596                ELSE
597                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))
598                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))
599                ENDIF
600             ENDIF   ! end if ( lconv )
601             !
602          END DO  ! end of ji loop
603       END DO     ! end of jj loop
604
605       !
606       ! calculate non-gradient components of the flux-gradient relationships
607       !
608! Stokes term in scalar flux, flux-gradient relationship
609       WHERE ( lconv )
610          zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln)
611          !
612          zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
613       ELSEWHERE
614          zsc_wth_1 = 2.0 * zwthav
615          !
616          zsc_ws_1 = 2.0 * zwsav
617       ENDWHERE
618
619
620       DO jj = 2, jpjm1
621          DO ji = 2, jpim1
622            IF ( lconv(ji,jj) ) THEN
623              DO jk = 2, imld(ji,jj)
624                 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
625                 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)
626                 !
627                 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)
628              END DO ! end jk loop
629            ELSE     ! else for if (lconv)
630 ! Stable conditions
631               DO jk = 2, ibld(ji,jj)
632                  zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
633                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
634                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
635                  !
636                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
637                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
638               END DO
639            ENDIF               ! endif for check on lconv
640
641          END DO  ! end of ji loop
642       END DO     ! end of jj loop
643
644! 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)
645       WHERE ( lconv )
646          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 )
647          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 )
648          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 )
649       ELSEWHERE
650          zsc_uw_1 = zustar**2
651          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln)
652       ENDWHERE
653       IF(ln_dia_osm) THEN
654          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
655          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
656       END IF
657       DO jj = 2, jpjm1
658          DO ji = 2, jpim1
659             IF ( lconv(ji,jj) ) THEN
660                DO jk = 2, imld(ji,jj)
661                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
662                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   &
663                        &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) &
664                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) )
665!
666                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       &
667                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
668                END DO   ! end jk loop
669             ELSE
670! Stable conditions
671                DO jk = 2, ibld(ji,jj) ! corrected to ibld
672                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
673                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       &
674                        &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
675                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
676                END DO   ! end jk loop
677             ENDIF
678          END DO  ! ji loop
679       END DO     ! jj loo
680
681! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
682
683       WHERE ( lconv )
684          zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
685          zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
686       ELSEWHERE
687          zsc_wth_1 = 0._wp
688          zsc_ws_1 = 0._wp
689       ENDWHERE
690
691       DO jj = 2, jpjm1
692          DO ji = 2, jpim1
693             IF (lconv(ji,jj) ) THEN
694                DO jk = 2, imld(ji,jj)
695                   zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
696                   ! calculate turbulent length scale
697                   zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
698                        &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) )
699                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
700                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
701                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
702                   ! non-gradient buoyancy terms
703                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
704                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
705                END DO
706             ELSE
707                DO jk = 2, ibld(ji,jj)
708                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
709                   ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
710                END DO
711             ENDIF
712          END DO   ! ji loop
713       END DO      ! jj loop
714
715       WHERE ( lconv )
716          zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
717          zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
718          zsc_vw_1 = 0._wp
719       ELSEWHERE
720         zsc_uw_1 = 0._wp
721         zsc_vw_1 = 0._wp
722       ENDWHERE
723
724       DO jj = 2, jpjm1
725          DO ji = 2, jpim1
726             IF ( lconv(ji,jj) ) THEN
727                DO jk = 2 , imld(ji,jj)
728                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
729                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     &
730                        &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   &
731                        &                                          * zsc_uw_2(ji,jj)                                    )
732                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
733                END DO  ! jk loop
734             ELSE
735             ! stable conditions
736                DO jk = 2, ibld(ji,jj)
737                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
738                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
739                END DO
740             ENDIF
741          END DO        ! ji loop
742       END DO           ! jj loop
743
744       IF(ln_dia_osm) THEN
745          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
746          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
747       END IF
748! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
749
750       WHERE ( lconv )
751          zsc_wth_1 = zwth0
752          zsc_ws_1 = zws0
753       ELSEWHERE
754          zsc_wth_1 = 2.0 * zwthav
755          zsc_ws_1 = zws0
756       ENDWHERE
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_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
763                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                &
764                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
765                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
766                       &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) )
767                  !
768                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  &
769                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
770                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
771                       &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) )
772               END DO
773            ELSE
774               DO jk = 2, ibld(ji,jj)
775                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
776                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
777                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
778               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
779                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
780               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
781               END DO
782            ENDIF
783          ENDDO    ! ji loop
784       END DO      ! jj loop
785
786       WHERE ( lconv )
787          zsc_uw_1 = zustar**2
788          zsc_vw_1 = ff_t * zustke * zhml
789       ELSEWHERE
790          zsc_uw_1 = zustar**2
791          zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
792          zsc_vw_1 = ff_t * zustke * zhbl
793          zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1
794       ENDWHERE
795
796       DO jj = 2, jpjm1
797          DO ji = 2, jpim1
798             IF ( lconv(ji,jj) ) THEN
799               DO jk = 2, imld(ji,jj)
800                  zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
801                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
802                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
803                       & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)
804                  !
805                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
806                       & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
807               END DO
808             ELSE
809               DO jk = 2, ibld(ji,jj)
810                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
811                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
812                  IF ( zznd_d <= 2.0 ) THEN
813                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
814                          &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
815                     !
816                  ELSE
817                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
818                          & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
819                     !
820                  ENDIF
821
822                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
823                       & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj)
824                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
825                       & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
826               END DO
827             ENDIF
828          END DO
829       END DO
830
831       IF(ln_dia_osm) THEN
832          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
833          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
834          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
835          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
836          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
837          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
838       END IF
839!
840! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
841
842      DO jj = 2, jpjm1
843         DO ji = 2, jpim1
844            IF ( lconv(ji,jj) ) THEN
845               DO jk = 2, ibld(ji,jj)
846                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
847                  IF ( znd >= 0.0 ) THEN
848                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
849                     ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
850                  ELSE
851                     ghamu(ji,jj,jk) = 0._wp
852                     ghamv(ji,jj,jk) = 0._wp
853                  ENDIF
854               END DO
855            ELSE
856               DO jk = 2, ibld(ji,jj)
857                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
858                  IF ( znd >= 0.0 ) THEN
859                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
860                     ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
861                  ELSE
862                     ghamu(ji,jj,jk) = 0._wp
863                     ghamv(ji,jj,jk) = 0._wp
864                  ENDIF
865               END DO
866            ENDIF
867         END DO
868      END DO
869
870       IF(ln_dia_osm) THEN
871          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
872          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
873       END IF
874      ! pynocline contributions
875       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small
876       zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln )
877       DO jj = 2, jpjm1
878          DO ji = 2, jpim1
879             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
880                DO jk= 2, ibld(ji,jj)
881                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
882                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
883                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
884                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
885                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk)
886                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
887                END DO
888             END IF
889          END DO
890       END DO
891
892! Entrainment contribution.
893
894       DO jj=2, jpjm1
895          DO ji = 2, jpim1
896            IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN
897               DO jk = 1, imld(ji,jj) - 1
898                  znd=gdepw_n(ji,jj,jk) / zhml(ji,jj)
899                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
900                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
901                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd
902                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd
903               END DO
904               DO jk = imld(ji,jj), ibld(ji,jj)
905                  znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
906                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
907                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
908                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd )
909                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd )
910                END DO
911             ENDIF
912
913             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
914             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
915             ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
916             ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
917          END DO       ! ji loop
918       END DO          ! jj loop
919
920       IF(ln_dia_osm) THEN
921          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
922          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
923          IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse )
924          IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse )
925          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
926          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
927          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
928       END IF
929       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
930       ! Need to put in code for contributions that are applied explicitly to
931       ! the prognostic variables
932       !  1. Entrainment flux
933       !
934       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
935
936
937
938       ! rotate non-gradient velocity terms back to model reference frame
939
940       DO jj = 2, jpjm1
941          DO ji = 2, jpim1
942             DO jk = 2, ibld(ji,jj)
943                ztemp = ghamu(ji,jj,jk)
944                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
945                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
946             END DO
947          END DO
948       END DO
949
950       IF(ln_dia_osm) THEN
951          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
952          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
953          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
954       END IF
955
956! KPP-style Ri# mixing
957       IF( ln_kpprimix) THEN
958          DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
959             DO jj = 1, jpjm1
960                DO ji = 1, jpim1   ! vector opt.
961                   z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
962                        &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
963                        &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
964                   z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
965                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
966                        &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
967                END DO
968             END DO
969          END DO
970      !
971         DO jk = 2, jpkm1
972            DO jj = 2, jpjm1
973               DO ji = 2, jpim1   ! vector opt.
974                  !                                          ! shear prod. at w-point weightened by mask
975                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
976                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
977                  !                                          ! local Richardson number
978                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
979                  zfri =  MIN( zri / rn_riinfty , 1.0_wp )
980                  zfri  = ( 1.0_wp - zfri * zfri )
981                  zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk)
982                END DO
983             END DO
984          END DO
985
986          DO jj = 2, jpjm1
987             DO ji = 2, jpim1
988                DO jk = ibld(ji,jj) + 1, jpkm1
989                   zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
990                   zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
991                END DO
992             END DO
993          END DO
994
995       END IF ! ln_kpprimix = .true.
996
997! KPP-style set diffusivity large if unstable below BL
998       IF( ln_convmix) THEN
999          DO jj = 2, jpjm1
1000             DO ji = 2, jpim1
1001                DO jk = ibld(ji,jj) + 1, jpkm1
1002                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
1003                END DO
1004             END DO
1005          END DO
1006       END IF ! ln_convmix = .true.
1007
1008       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
1009       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
1010
1011       ! GN 25/8: need to change tmask --> wmask
1012
1013     DO jk = 2, jpkm1
1014         DO jj = 2, jpjm1
1015             DO ji = 2, jpim1
1016                p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1017                p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
1018             END DO
1019         END DO
1020     END DO
1021      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid  (sign unchanged), needed to caclulate gham[uv] on u and v grids
1022     CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1.,   &
1023      &                  ghamu, 'W', 1. , ghamv, 'W', 1. )
1024       DO jk = 2, jpkm1
1025           DO jj = 2, jpjm1
1026               DO ji = 2, jpim1
1027                  ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
1028                     &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1029
1030                  ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
1031                      &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1032
1033                  ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk)
1034                  ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk)
1035               END DO
1036           END DO
1037        END DO
1038        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged)
1039        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged)
1040        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1.,   &
1041         &                  ghamu, 'U', -1. , ghamv, 'V', -1. )
1042
1043      IF(ln_dia_osm) THEN
1044         SELECT CASE (nn_osm_wave)
1045         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1).
1046         CASE(0:1)
1047            IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift
1048            IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift
1049            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
1050         ! Stokes drift read in from sbcwave  (=2).
1051         CASE(2)
1052            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift
1053            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift
1054            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period
1055            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height
1056            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
1057            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum
1058            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10
1059            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* &
1060                 & SQRT(ut0sd**2 + vt0sd**2 ) )
1061         END SELECT
1062         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )            ! <Tw_NL>
1063         IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams )            ! <Sw_NL>
1064         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )            ! <uw_NL>
1065         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )            ! <vw_NL>
1066         IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 )            ! <Tw_0>
1067         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0>
1068         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth
1069         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k
1070         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base
1071         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base
1072         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base
1073         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base
1074         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base
1075         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth
1076         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth
1077         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth
1078         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points
1079         IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc )         ! convective velocity scale
1080         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale
1081         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale
1082         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale
1083         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir #
1084         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
1085         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
1086         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine
1087         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine
1088         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine
1089         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine
1090         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine
1091         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux
1092         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux
1093         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux
1094         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux
1095         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML
1096      END IF
1097
1098CONTAINS
1099
1100
1101   ! Alan: do we need zb?
1102   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv )
1103     !!---------------------------------------------------------------------
1104     !!                   ***  ROUTINE zdf_vertical_average  ***
1105     !!
1106     !! ** Purpose : Determines vertical averages from surface to jnlev.
1107     !!
1108     !! ** Method  : Averages are calculated from the surface to jnlev.
1109     !!              The external level used to calculate differences is ibld+ibld_ext
1110     !!
1111     !!----------------------------------------------------------------------
1112
1113        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over.
1114
1115        ! Alan: do we need zb?
1116        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity
1117        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components
1118        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
1119        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components.
1120
1121        INTEGER :: jk, ji, jj
1122        REAL(wp) :: zthick, zthermal, zbeta
1123
1124
1125        zt   = 0._wp
1126        zs   = 0._wp
1127        zu   = 0._wp
1128        zv   = 0._wp
1129        DO jj = 2, jpjm1                                 !  Vertical slab
1130         DO ji = 2, jpim1
1131            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1132            zbeta    = rab_n(ji,jj,1,jp_sal)
1133               ! average over depth of boundary layer
1134            zthick = epsln
1135            DO jk = 2, jnlev_av(ji,jj)
1136               zthick = zthick + e3t_n(ji,jj,jk)
1137               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
1138               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
1139               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) &
1140                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
1141                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
1142               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) &
1143                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
1144                     &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
1145            END DO
1146            zt(ji,jj) = zt(ji,jj) / zthick
1147            zs(ji,jj) = zs(ji,jj) / zthick
1148            zu(ji,jj) = zu(ji,jj) / zthick
1149            zv(ji,jj) = zv(ji,jj) / zthick
1150            ! Alan: do we need zb?
1151            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem)
1152            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal)
1153            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) &
1154                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) )
1155            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) &
1156                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) )
1157            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
1158         END DO
1159        END DO
1160   END SUBROUTINE zdf_osm_vertical_average
1161
1162   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
1163     !!---------------------------------------------------------------------
1164     !!                   ***  ROUTINE zdf_velocity_rotation  ***
1165     !!
1166     !! ** Purpose : Rotates frame of reference of averaged velocity components.
1167     !!
1168     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w
1169     !!
1170     !!----------------------------------------------------------------------
1171
1172        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle
1173        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current
1174        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline
1175
1176        INTEGER :: ji, jj
1177        REAL(wp) :: ztemp
1178
1179        DO jj = 2, jpjm1
1180           DO ji = 2, jpim1
1181              ztemp = zu(ji,jj)
1182              zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
1183              zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
1184              ztemp = zdu(ji,jj)
1185              zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
1186              zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj)
1187           END DO
1188        END DO
1189    END SUBROUTINE zdf_osm_velocity_rotation
1190
1191    SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz )
1192     !!---------------------------------------------------------------------
1193     !!                   ***  ROUTINE zdf_osm_external_gradients  ***
1194     !!
1195     !! ** Purpose : Calculates the gradients below the OSBL
1196     !!
1197     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient.
1198     !!
1199     !!----------------------------------------------------------------------
1200
1201     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy.
1202
1203     INTEGER :: jj, ji, jkb, jkb1
1204     REAL(wp) :: zthermal, zbeta
1205
1206
1207     DO jj = 2, jpjm1
1208        DO ji = 2, jpim1
1209           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1210              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1211              zbeta    = rab_n(ji,jj,1,jp_sal)
1212              jkb = ibld(ji,jj) + ibld_ext
1213              jkb1 = MIN(jkb + 1, mbkt(ji,jj))
1214              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
1215                   &   / e3t_n(ji,jj,ibld(ji,jj))
1216              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
1217                   &   / e3t_n(ji,jj,ibld(ji,jj))
1218              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
1219           END IF
1220        END DO
1221     END DO
1222    END SUBROUTINE zdf_osm_external_gradients
1223
1224    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz )
1225
1226     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline
1227
1228     INTEGER :: jk, jj, ji
1229     REAL(wp) :: ztgrad, zsgrad, zbgrad
1230     REAL(wp) :: zgamma_b_nd, zgamma_c, znd
1231     REAL(wp) :: zzeta_s=0.3
1232
1233     DO jj = 2, jpjm1
1234        DO ji = 2, jpim1
1235           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1236              IF ( lconv(ji,jj) ) THEN  ! convective conditions
1237                 IF ( zdb_bl(ji,jj) > 0._wp .AND. zdbdz_ext(ji,jj) > 0._wp ) THEN
1238                    ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj)
1239                    zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj)
1240                    zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj)
1241                    zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1242                    zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
1243                         & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
1244                    DO jk = 2, ibld(ji,jj)+ibld_ext
1245                       znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1246                       IF ( znd <= zzeta_s ) THEN
1247                          zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * &
1248                               & EXP( -6.0 * ( znd -zzeta_s )**2 )
1249                          zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * &
1250                               & EXP( -6.0 * ( znd -zzeta_s )**2 )
1251                          zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * &
1252                               & EXP( -6.0 * ( znd -zzeta_s )**2 )
1253                       ELSE
1254                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1255                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1256                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1257                       ENDIF
1258                    END DO
1259                 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing
1260              ELSE
1261                 ! stable conditions
1262                 ! if pycnocline profile only defined when depth steady of increasing.
1263                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady.
1264                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1265                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
1266                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
1267                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
1268                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
1269                          DO jk = 2, ibld(ji,jj)
1270                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1271                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1272                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1273                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1274                          END DO
1275                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
1276                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
1277                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
1278                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
1279                          DO jk = 2, ibld(ji,jj)
1280                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
1281                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1282                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1283                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1284                          END DO
1285                       ENDIF ! IF (zhol >=0.5)
1286                    ENDIF    ! IF (zdb_bl> 0.)
1287                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
1288              ENDIF          ! IF (lconv)
1289           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1290        END DO
1291     END DO
1292
1293    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
1294
1295    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
1296      !!---------------------------------------------------------------------
1297      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  ***
1298      !!
1299      !! ** Purpose : Calculates velocity shear in the pycnocline
1300      !!
1301      !! ** Method  :
1302      !!
1303      !!----------------------------------------------------------------------
1304
1305      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
1306
1307      INTEGER :: jk, jj, ji
1308      REAL(wp) :: zugrad, zvgrad, znd
1309      REAL(wp) :: zzeta_v = 0.45
1310      !
1311      DO jj = 2, jpjm1
1312         DO ji = 2, jpim1
1313            !
1314            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1315               IF ( lconv (ji,jj) ) THEN
1316                  ! Unstable conditions
1317                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
1318                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
1319                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
1320                  !Alan is this right?
1321                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
1322                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
1323                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
1324                       &      )/ (zdh(ji,jj)  + epsln )
1325                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
1326                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
1327                     IF ( znd <= 0.0 ) THEN
1328                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
1329                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
1330                     ELSE
1331                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
1332                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
1333                     ENDIF
1334                  END DO
1335               ELSE
1336                  ! stable conditions
1337                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
1338                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
1339                  DO jk = 2, ibld(ji,jj)
1340                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1341                     IF ( znd < 1.0 ) THEN
1342                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
1343                     ELSE
1344                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
1345                     ENDIF
1346                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
1347                  END DO
1348               ENDIF
1349               !
1350            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1351         END DO
1352      END DO
1353    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
1354
1355    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
1356     !!---------------------------------------------------------------------
1357     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
1358     !!
1359     !! ** Purpose : Calculates the rate at which hbl changes.
1360     !!
1361     !! ** Method  :
1362     !!
1363     !!----------------------------------------------------------------------
1364
1365    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl
1366
1367    INTEGER :: jj, ji
1368    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert
1369    REAL(wp) :: zvel_max, zwb_min
1370    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
1371    REAL(wp) :: zzeta_m = 0.3
1372    REAL(wp) :: zgamma_c = 2.0
1373    REAL(wp) :: zdhoh = 0.1
1374    REAL(wp) :: alpha_bc = 0.5
1375
1376    DO jj = 2, jpjm1
1377       DO ji = 2, jpim1
1378          IF ( lconv(ji,jj) ) THEN    ! Convective
1379             ! Alan is this right?  Yes, it's a bit different from the previous relationship
1380             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) &
1381             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1382             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
1383             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
1384             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
1385             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
1386             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
1387                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
1388
1389             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
1390                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
1391                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
1392                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1393             !
1394             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj)
1395                  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) / &
1396                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1397                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1398                ! added ajgn 23 July as temporay fix
1399             zdhdt_2(ji,jj) = 0._wp
1400
1401                ! commented out ajgn 23 July as temporay fix
1402!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1403! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary.
1404!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1405!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj)
1406!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1407!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj)
1408!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min )
1409!                     ! Alan no idea what this should be?
1410!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) &
1411!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) &
1412!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj)
1413!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
1414!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN
1415!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj)
1416!                 ENDIF
1417            ELSE                        ! Stable
1418                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
1419                zdhdt_2(ji,jj) = 0._wp
1420                IF ( zdhdt(ji,jj) < 0._wp ) THEN
1421                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
1422                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
1423                ELSE
1424                    zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
1425                ENDIF
1426                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert
1427            ENDIF
1428         END DO
1429      END DO
1430    END SUBROUTINE zdf_osm_calculate_dhdt
1431
1432    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
1433     !!---------------------------------------------------------------------
1434     !!                   ***  ROUTINE zdf_osm_timestep_hbl  ***
1435     !!
1436     !! ** Purpose : Increments hbl.
1437     !!
1438     !! ** Method  : If thechange in hbl exceeds one model level the change is
1439     !!              is calculated by moving down the grid, changing the buoyancy
1440     !!              jump. This is to ensure that the change in hbl does not
1441     !!              overshoot a stable layer.
1442     !!
1443     !!----------------------------------------------------------------------
1444
1445
1446    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl.
1447
1448    INTEGER :: jk, jj, ji, jm
1449    REAL(wp) :: zhbl_s, zvel_max, zdb
1450    REAL(wp) :: zthermal, zbeta
1451
1452     DO jj = 2, jpjm1
1453         DO ji = 2, jpim1
1454           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
1455!
1456! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
1457!
1458              zhbl_s = hbl(ji,jj)
1459              jm = imld(ji,jj)
1460              zthermal = rab_n(ji,jj,1,jp_tem)
1461              zbeta = rab_n(ji,jj,1,jp_sal)
1462
1463
1464              IF ( lconv(ji,jj) ) THEN
1465!unstable
1466                    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) / &
1467                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1468                 DO jk = imld(ji,jj), ibld(ji,jj)
1469                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
1470                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
1471                         &  0.0 ) + zvel_max
1472
1473                      zhbl_s = zhbl_s + MIN( &
1474                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
1475                        & e3w_n(ji,jj,jm) )
1476                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1477
1478                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
1479                 END DO
1480                 hbl(ji,jj) = zhbl_s
1481                 ibld(ji,jj) = jm
1482             ELSE
1483! stable
1484                 DO jk = imld(ji,jj), ibld(ji,jj)
1485                    zdb = MAX( &
1486                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&
1487                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
1488                         & 0.0 ) + &
1489             &       2.0 * zvstr(ji,jj)**2 / zhbl_s
1490
1491                    ! Alan is thuis right? I have simply changed hbli to hbl
1492                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
1493                    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) ) ) * &
1494               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
1495                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
1496                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
1497
1498                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1499                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
1500                 END DO
1501             ENDIF   ! IF ( lconv )
1502             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
1503             ibld(ji,jj) = MAX(jm, 4 )
1504           ELSE
1505! change zero or one model level.
1506             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
1507           ENDIF
1508           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
1509         END DO
1510      END DO
1511
1512    END SUBROUTINE zdf_osm_timestep_hbl
1513
1514    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
1515      !!---------------------------------------------------------------------
1516      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  ***
1517      !!
1518      !! ** Purpose : Calculates thickness of the pycnocline
1519      !!
1520      !! ** Method  : The thickness is calculated from a prognostic equation
1521      !!              that relaxes the pycnocine thickness to a diagnostic
1522      !!              value. The time change is calculated assuming the
1523      !!              thickness relaxes exponentially. This is done to deal
1524      !!              with large timesteps.
1525      !!
1526      !!----------------------------------------------------------------------
1527
1528      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness.
1529      !
1530      INTEGER :: jj, ji
1531      INTEGER :: inhml
1532      REAL(wp) :: zari, ztau, zddhdt
1533
1534
1535      DO jj = 2, jpjm1
1536         DO ji = 2, jpim1
1537            IF ( lconv(ji,jj) ) THEN
1538               IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1539                  IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
1540                     zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1541                          (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1542                  ELSE                                                     ! unstable
1543                     zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1544                          (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1545                  ENDIF
1546                  ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
1547                  zddhdt = zari * hbl(ji,jj)
1548               ELSE
1549                  ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1550                  zddhdt = 0.2 * hbl(ji,jj)
1551               ENDIF
1552               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1553               ! Alan: this hml is never defined or used
1554            ELSE   ! IF (lconv)
1555               ztau = hbl(ji,jj) / zvstr(ji,jj)
1556               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
1557                  ! boundary layer deepening
1558                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1559                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
1560                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
1561                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
1562                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj)
1563                  ELSE
1564                     zddhdt = 0.2 * hbl(ji,jj)
1565                  ENDIF
1566               ELSE     ! IF(dhdt < 0)
1567                  zddhdt = 0.2 * hbl(ji,jj)
1568               ENDIF    ! IF (dhdt >= 0)
1569               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1570               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt       ! can be a problem with dh>hbl for rapid collapse
1571               ! Alan: this hml is never defined or used -- do we need it?
1572            ENDIF       ! IF (lconv)
1573
1574            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
1575            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
1576            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
1577            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
1578            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1579         END DO
1580      END DO
1581
1582    END SUBROUTINE zdf_osm_pycnocline_thickness
1583
1584END SUBROUTINE zdf_osm
1585
1586
1587   SUBROUTINE zdf_osm_init
1588     !!----------------------------------------------------------------------
1589     !!                  ***  ROUTINE zdf_osm_init  ***
1590     !!
1591     !! ** Purpose :   Initialization of the vertical eddy diffivity and
1592     !!      viscosity when using a osm turbulent closure scheme
1593     !!
1594     !! ** Method  :   Read the namosm namelist and check the parameters
1595     !!      called at the first timestep (nit000)
1596     !!
1597     !! ** input   :   Namlist namosm
1598     !!----------------------------------------------------------------------
1599     INTEGER  ::   ios            ! local integer
1600     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1601     REAL z1_t2
1602     !!
1603     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
1604          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
1605          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave
1606     !!----------------------------------------------------------------------
1607     !
1608     REWIND( numnam_ref )              ! Namelist namzdf_osm in reference namelist : Osmosis ML model
1609     READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
1610901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist', lwp )
1611
1612     REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
1613     READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
1614902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist', lwp )
1615     IF(lwm) WRITE ( numond, namzdf_osm )
1616
1617     IF(lwp) THEN                    ! Control print
1618        WRITE(numout,*)
1619        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
1620        WRITE(numout,*) '~~~~~~~~~~~~'
1621        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters'
1622        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la
1623        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
1624        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
1625        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes
1626        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
1627        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave
1628        SELECT CASE (nn_osm_wave)
1629        CASE(0)
1630           WRITE(numout,*) '     calculated assuming constant La#=0.3'
1631        CASE(1)
1632           WRITE(numout,*) '     calculated from Pierson Moskowitz wind-waves'
1633        CASE(2)
1634           WRITE(numout,*) '     calculated from ECMWF wave fields'
1635        END SELECT
1636        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
1637        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix
1638        WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
1639        WRITE(numout,*) '     maximum shear diffusivity at Rig = 0    (m2/s) rn_difri = ', rn_difri
1640        WRITE(numout,*) '     Use large mixing below BL when unstable       ln_convmix = ', ln_convmix
1641        WRITE(numout,*) '     diffusivity when unstable below BL     (m2/s) rn_difconv = ', rn_difconv
1642     ENDIF
1643
1644     !                              ! allocate zdfosm arrays
1645     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
1646
1647
1648     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle
1649
1650
1651     IF( ln_zdfddm) THEN
1652        IF(lwp) THEN
1653           WRITE(numout,*)
1654           WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
1655           WRITE(numout,*) '    CAUTION : done in routine zdfosm, not in routine zdfddm '
1656        ENDIF
1657     ENDIF
1658
1659
1660     !set constants not in namelist
1661     !-----------------------------
1662
1663     IF(lwp) THEN
1664        WRITE(numout,*)
1665     ENDIF
1666
1667     IF (nn_osm_wave == 0) THEN
1668        dstokes(:,:) = rn_osm_dstokes
1669     END IF
1670
1671     ! Horizontal average : initialization of weighting arrays
1672     ! -------------------
1673
1674     SELECT CASE ( nn_ave )
1675
1676     CASE ( 0 )                ! no horizontal average
1677        IF(lwp) WRITE(numout,*) '          no horizontal average on avt'
1678        IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
1679        ! weighting mean arrays etmean
1680        !           ( 1  1 )
1681        ! avt = 1/4 ( 1  1 )
1682        !
1683        etmean(:,:,:) = 0.e0
1684
1685        DO jk = 1, jpkm1
1686           DO jj = 2, jpjm1
1687              DO ji = 2, jpim1   ! vector opt.
1688                 etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
1689                      &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
1690                      &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
1691              END DO
1692           END DO
1693        END DO
1694
1695     CASE ( 1 )                ! horizontal average
1696        IF(lwp) WRITE(numout,*) '          horizontal average on avt'
1697        ! weighting mean arrays etmean
1698        !           ( 1/2  1  1/2 )
1699        ! avt = 1/8 ( 1    2  1   )
1700        !           ( 1/2  1  1/2 )
1701        etmean(:,:,:) = 0.e0
1702
1703        DO jk = 1, jpkm1
1704           DO jj = 2, jpjm1
1705              DO ji = 2, jpim1   ! vector opt.
1706                 etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
1707                      & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
1708                      &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
1709                      &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
1710                      &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
1711                      &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
1712              END DO
1713           END DO
1714        END DO
1715
1716     CASE DEFAULT
1717        WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
1718        CALL ctl_stop( ctmp1 )
1719
1720     END SELECT
1721
1722     ! Initialization of vertical eddy coef. to the background value
1723     ! -------------------------------------------------------------
1724     DO jk = 1, jpk
1725        avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1726     END DO
1727
1728     ! zero the surface flux for non local term and osm mixed layer depth
1729     ! ------------------------------------------------------------------
1730     ghamt(:,:,:) = 0.
1731     ghams(:,:,:) = 0.
1732     ghamu(:,:,:) = 0.
1733     ghamv(:,:,:) = 0.
1734     !
1735     IF( lwxios ) THEN
1736        CALL iom_set_rstw_var_active('wn')
1737        CALL iom_set_rstw_var_active('hbl')
1738        CALL iom_set_rstw_var_active('hbli')
1739     ENDIF
1740   END SUBROUTINE zdf_osm_init
1741
1742
1743   SUBROUTINE osm_rst( kt, cdrw )
1744     !!---------------------------------------------------------------------
1745     !!                   ***  ROUTINE osm_rst  ***
1746     !!
1747     !! ** Purpose :   Read or write BL fields in restart file
1748     !!
1749     !! ** Method  :   use of IOM library. If the restart does not contain
1750     !!                required fields, they are recomputed from stratification
1751     !!----------------------------------------------------------------------
1752
1753     INTEGER, INTENT(in) :: kt
1754     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1755
1756     INTEGER ::   id1, id2   ! iom enquiry index
1757     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1758     INTEGER  ::   iiki, ikt ! local integer
1759     REAL(wp) ::   zhbf           ! tempory scalars
1760     REAL(wp) ::   zN2_c           ! local scalar
1761     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
1762     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
1763     !!----------------------------------------------------------------------
1764     !
1765     !!-----------------------------------------------------------------------------
1766     ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
1767     !!-----------------------------------------------------------------------------
1768     IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
1769        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. )
1770        IF( id1 > 0 ) THEN                       ! 'wn' exists; read
1771           CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios )
1772           WRITE(numout,*) ' ===>>>> :  wn read from restart file'
1773        ELSE
1774           wn(:,:,:) = 0._wp
1775           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
1776        END IF
1777
1778        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
1779        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. )
1780        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return
1781           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
1782           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  )
1783           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file'
1784           RETURN
1785        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate
1786           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
1787        END IF
1788     END IF
1789
1790     !!-----------------------------------------------------------------------------
1791     ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
1792     !!-----------------------------------------------------------------------------
1793     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
1794        IF(lwp) WRITE(numout,*) '---- osm-rst ----'
1795         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios )
1796         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios )
1797         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh, ldxios = lwxios )
1798        RETURN
1799     END IF
1800
1801     !!-----------------------------------------------------------------------------
1802     ! Getting hbl, no restart file with hbl, so calculate from surface stratification
1803     !!-----------------------------------------------------------------------------
1804     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
1805     ! w-level of the mixing and mixed layers
1806     CALL eos_rab( tsn, rab_n )
1807     CALL bn2(tsn, rab_n, rn2)
1808     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
1809     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
1810     zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
1811     !
1812     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
1813     DO jk = 1, jpkm1
1814        DO jj = 1, jpj              ! Mixed layer level: w-level
1815           DO ji = 1, jpi
1816              ikt = mbkt(ji,jj)
1817              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
1818              IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
1819           END DO
1820        END DO
1821     END DO
1822     !
1823     DO jj = 1, jpj
1824        DO ji = 1, jpi
1825           iiki = MAX(4,imld_rst(ji,jj))
1826           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth
1827           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth
1828        END DO
1829     END DO
1830     WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
1831     wn(:,:,:) = 0._wp
1832     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
1833   END SUBROUTINE osm_rst
1834
1835
1836   SUBROUTINE tra_osm( kt )
1837      !!----------------------------------------------------------------------
1838      !!                  ***  ROUTINE tra_osm  ***
1839      !!
1840      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
1841      !!
1842      !! ** Method  :   ???
1843      !!----------------------------------------------------------------------
1844      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
1845      !!----------------------------------------------------------------------
1846      INTEGER, INTENT(in) :: kt
1847      INTEGER :: ji, jj, jk
1848      !
1849      IF( kt == nit000 ) THEN
1850         IF(lwp) WRITE(numout,*)
1851         IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
1852         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1853      ENDIF
1854
1855      IF( l_trdtra )   THEN                    !* Save ta and sa trends
1856         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
1857         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
1858      ENDIF
1859
1860      DO jk = 1, jpkm1
1861         DO jj = 2, jpjm1
1862            DO ji = 2, jpim1
1863               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
1864                  &                 - (  ghamt(ji,jj,jk  )  &
1865                  &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
1866               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
1867                  &                 - (  ghams(ji,jj,jk  )  &
1868                  &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
1869            END DO
1870         END DO
1871      END DO
1872
1873      ! save the non-local tracer flux trends for diagnostics
1874      IF( l_trdtra )   THEN
1875         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
1876         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
1877
1878         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt )
1879         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds )
1880         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
1881      ENDIF
1882
1883      IF(ln_ctl) THEN
1884         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   &
1885         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
1886      ENDIF
1887      !
1888   END SUBROUTINE tra_osm
1889
1890
1891   SUBROUTINE trc_osm( kt )          ! Dummy routine
1892      !!----------------------------------------------------------------------
1893      !!                  ***  ROUTINE trc_osm  ***
1894      !!
1895      !! ** Purpose :   compute and add to the passive tracer trend the non-local
1896      !!                 passive tracer flux
1897      !!
1898      !!
1899      !! ** Method  :   ???
1900      !!----------------------------------------------------------------------
1901      !
1902      !!----------------------------------------------------------------------
1903      INTEGER, INTENT(in) :: kt
1904      WRITE(*,*) 'trc_osm: Not written yet', kt
1905   END SUBROUTINE trc_osm
1906
1907
1908   SUBROUTINE dyn_osm( kt )
1909      !!----------------------------------------------------------------------
1910      !!                  ***  ROUTINE dyn_osm  ***
1911      !!
1912      !! ** Purpose :   compute and add to the velocity trend the non-local flux
1913      !! copied/modified from tra_osm
1914      !!
1915      !! ** Method  :   ???
1916      !!----------------------------------------------------------------------
1917      INTEGER, INTENT(in) ::   kt   !
1918      !
1919      INTEGER :: ji, jj, jk   ! dummy loop indices
1920      !!----------------------------------------------------------------------
1921      !
1922      IF( kt == nit000 ) THEN
1923         IF(lwp) WRITE(numout,*)
1924         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
1925         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1926      ENDIF
1927      !code saving tracer trends removed, replace with trdmxl_oce
1928
1929      DO jk = 1, jpkm1           ! add non-local u and v fluxes
1930         DO jj = 2, jpjm1
1931            DO ji = 2, jpim1
1932               ua(ji,jj,jk) =  ua(ji,jj,jk)                      &
1933                  &                 - (  ghamu(ji,jj,jk  )  &
1934                  &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
1935               va(ji,jj,jk) =  va(ji,jj,jk)                      &
1936                  &                 - (  ghamv(ji,jj,jk  )  &
1937                  &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
1938            END DO
1939         END DO
1940      END DO
1941      !
1942      ! code for saving tracer trends removed
1943      !
1944   END SUBROUTINE dyn_osm
1945
1946   !!======================================================================
1947
1948END MODULE zdfosm
Note: See TracBrowser for help on using the repository browser.