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

Last change on this file since 11145 was 11145, checked in by cguiavarch, 21 months ago

Merge changes from /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser:r11139-11143

File size: 101.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) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
370           END DO
371        END DO
372     ! Assume Pierson-Moskovitz wind-wave spectrum
373     CASE(1)
374        DO jj = 2, jpjm1
375           DO ji = 2, jpim1
376              ! Use wind speed wndm included in sbc_oce module
377              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
378              dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav
379           END DO
380        END DO
381     ! Use ECMWF wave fields as output from SBCWAVE
382     CASE(2)
383        zfac =  2.0_wp * rpi / 16.0_wp
384
385        DO jj = 2, jpjm1
386           DO ji = 2, jpim1
387              ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
388              !    The coefficient 0.8 gives La=0.3  in this situation.
389              ! It could represent the effects of the spread of wave directions
390              ! around the mean wind. The effect of this adjustment needs to be tested.
391              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)
392              zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8)
393              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 !
394           END DO
395        END DO
396     END SELECT
397
398     ! Langmuir velocity scale (zwstrl), La # (zla)
399     ! mixed scale (zvstr), convective velocity scale (zwstrc)
400     DO jj = 2, jpjm1
401        DO ji = 2, jpim1
402           ! Langmuir velocity scale (zwstrl), at T-point
403           zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
404           zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)
405           IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))
406           ! Velocity scale that tends to zustar for large Langmuir numbers
407           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + &
408                & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
409
410           ! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
411           ! Note zustke and zwstrl are not amended.
412           !
413           ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
414           IF ( zwbav(ji,jj) > 0.0) THEN
415              zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
416              zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
417              lconv(ji,jj) = .TRUE.
418           ELSE
419              zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln )
420              lconv(ji,jj) = .FALSE.
421           ENDIF
422        END DO
423     END DO
424
425     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
426     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
427     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
428     ! BL must be always 4 levels deep.
429      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) )
430      ibld(:,:) = 4
431      DO jk = 5, jpkm1
432         DO jj = 2, jpjm1
433            DO ji = 2, jpim1
434               IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
435                  ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
436               ENDIF
437            END DO
438         END DO
439      END DO
440
441      DO jj = 2, jpjm1
442         DO ji = 2, jpim1
443            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
444            imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 ))
445            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
446         END DO
447      END DO
448      ! Averages over well-mixed and boundary layer
449      ! Alan: do we need zb_nl?, zb_ml?
450      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl)
451      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
452! External gradient
453      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )
454
455! Rate of change of hbl
456      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
457      ! Calculate averages over depth of boundary layer
458      imld = ibld           ! use imld to hold previous blayer index
459      ibld(:,:) = 4
460
461         DO jj = 2, jpjm1
462            DO ji = 2, jpim1
463               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
464               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:))
465               !zdhdt(ji,jj) = MIN(zdhdt(ji,jj), (zhbl_t(ji,jj) - hbl(ji,jj))/rn_rdt + wn(ji,jj,ibld(ji,jj)))
466            END DO
467         END DO
468      ! adjustment to represent limiting by ocean bottom
469
470      DO jk = 4, jpkm1
471         DO jj = 2, jpjm1
472            DO ji = 2, jpim1
473               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
474                  ibld(ji,jj) = jk
475               ENDIF
476            END DO
477         END DO
478      END DO
479
480!
481! Step through model levels taking account of buoyancy change to determine the effect on dhdt
482!
483      CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
484      ! Alan: do we need zb_ml?
485      CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )
486!
487!
488      CALL zdf_osm_pycnocline_thickness( dh, zdh )
489      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms.
490!
491    ! Average over the depth of the mixed layer in the convective boundary layer
492    ! Alan: do we need zb_ml?
493    CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )
494    ! rotate mean currents and changes onto wind align co-ordinates
495    !
496    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
497    CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
498     zuw_bse = 0._wp
499     zvw_bse = 0._wp
500     zwth_ent = 0._wp
501     zws_ent = 0._wp
502     DO jj = 2, jpjm1
503        DO ji = 2, jpim1
504           IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN
505              IF ( lconv(ji,jj) ) THEN
506             zuw_bse(ji,jj) = -0.0075*((zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdu_ml(ji,jj) + &
507                      &                    1.5*zustar(ji,jj)**2*(zhbl(ji,jj)-zhml(ji,jj)) )/ &
508                      &                     ( zhml(ji,jj)*MIN(zla(ji,jj)**(8./3.),1.) + epsln)
509            zvw_bse(ji,jj) = 0.01*(-(zvstr(ji,jj)**3+0.5*zwstrc(ji,jj)**3)**pthird*zdv_ml(ji,jj)+ &
510                      &                    2.0*ff_t(ji,jj)*zustke(ji,jj)*dstokes(ji,jj)*zla(ji,jj))
511                 IF ( zdb_ml(ji,jj) > 0._wp ) THEN
512                    zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
513                    zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
514                 ENDIF
515              ELSE
516                 zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )
517                 zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )
518              ENDIF
519           ENDIF
520        END DO
521     END DO
522
523      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
524      !  Pycnocline gradients for scalars and velocity
525      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
526
527      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )
528      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc )
529      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )
530       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
531       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
532       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
533
534       DO jj = 2, jpjm1
535          DO ji = 2, jpim1
536             IF ( lconv(ji,jj) ) THEN
537               zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
538               zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
539               zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
540               zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
541               zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third
542               zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln )
543             ELSE
544               zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
545               zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
546             END IF
547          END DO
548       END DO
549!
550       DO jj = 2, jpjm1
551          DO ji = 2, jpim1
552             IF ( lconv(ji,jj) ) THEN
553                DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
554                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
555                    !
556                    zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5
557                    !
558                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) &
559                         &            *                                      ( 1.0 -               0.5 * zznd_ml**2 )
560                END DO
561                ! pycnocline - if present linear profile
562                IF ( zdh(ji,jj) > 0._wp ) THEN
563                   zgamma_b = 6.0
564                   DO jk = imld(ji,jj)+1 , ibld(ji,jj)
565                       zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
566                       !
567                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
568                       !
569                       zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
570                   END DO
571                   IF ( ibld_ext == 0 ) THEN
572                       zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
573                       zviscos(ji,jj,ibld(ji,jj)) = 0._wp
574                   ELSE
575                       zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )
576                       zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw_n(ji, jj, ibld(ji,jj)-1) )
577                   ENDIF
578                ENDIF
579                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out
580                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj))
581                zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj))
582                ! could be taken out, take account of entrainment represents as a diffusivity
583                ! should remove w from here, represents entrainment
584             ELSE
585             ! stable conditions
586                DO jk = 2, ibld(ji,jj)
587                   zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
588                   zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
589                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
590                END DO
591
592                IF ( ibld_ext == 0 ) THEN
593                   zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
594                   zviscos(ji,jj,ibld(ji,jj)) = 0._wp
595                ELSE
596                   zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))
597                   zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w_n(ji, jj, ibld(ji,jj))
598                ENDIF
599             ENDIF   ! end if ( lconv )
600             !
601          END DO  ! end of ji loop
602       END DO     ! end of jj loop
603
604       !
605       ! calculate non-gradient components of the flux-gradient relationships
606       !
607! Stokes term in scalar flux, flux-gradient relationship
608       WHERE ( lconv )
609          zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln)
610          !
611          zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
612       ELSEWHERE
613          zsc_wth_1 = 2.0 * zwthav
614          !
615          zsc_ws_1 = 2.0 * zwsav
616       ENDWHERE
617
618
619       DO jj = 2, jpjm1
620          DO ji = 2, jpim1
621            IF ( lconv(ji,jj) ) THEN
622              DO jk = 2, imld(ji,jj)
623                 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
624                 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)
625                 !
626                 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)
627              END DO ! end jk loop
628            ELSE     ! else for if (lconv)
629 ! Stable conditions
630               DO jk = 2, ibld(ji,jj)
631                  zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
632                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
633                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
634                  !
635                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
636                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
637               END DO
638            ENDIF               ! endif for check on lconv
639
640          END DO  ! end of ji loop
641       END DO     ! end of jj loop
642
643! 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)
644       WHERE ( lconv )
645          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 )
646          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 )
647          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 )
648       ELSEWHERE
649          zsc_uw_1 = zustar**2
650          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln)
651       ENDWHERE
652       IF(ln_dia_osm) THEN
653          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
654          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
655       END IF
656       DO jj = 2, jpjm1
657          DO ji = 2, jpim1
658             IF ( lconv(ji,jj) ) THEN
659                DO jk = 2, imld(ji,jj)
660                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
661                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   &
662                        &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) &
663                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) )
664!
665                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       &
666                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
667                END DO   ! end jk loop
668             ELSE
669! Stable conditions
670                DO jk = 2, ibld(ji,jj) ! corrected to ibld
671                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
672                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       &
673                        &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
674                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
675                END DO   ! end jk loop
676             ENDIF
677          END DO  ! ji loop
678       END DO     ! jj loo
679
680! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
681
682       WHERE ( lconv )
683          zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
684          zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
685       ELSEWHERE
686          zsc_wth_1 = 0._wp
687          zsc_ws_1 = 0._wp
688       ENDWHERE
689
690       DO jj = 2, jpjm1
691          DO ji = 2, jpim1
692             IF (lconv(ji,jj) ) THEN
693                DO jk = 2, imld(ji,jj)
694                   zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
695                   ! calculate turbulent length scale
696                   zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
697                        &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) )
698                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
699                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
700                   zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
701                   ! non-gradient buoyancy terms
702                   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 )
703                   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 )
704                END DO
705             ELSE
706                DO jk = 2, ibld(ji,jj)
707                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
708                   ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
709                END DO
710             ENDIF
711          END DO   ! ji loop
712       END DO      ! jj loop
713
714       WHERE ( lconv )
715          zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
716          zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
717          zsc_vw_1 = 0._wp
718       ELSEWHERE
719         zsc_uw_1 = 0._wp
720         zsc_vw_1 = 0._wp
721       ENDWHERE
722
723       DO jj = 2, jpjm1
724          DO ji = 2, jpim1
725             IF ( lconv(ji,jj) ) THEN
726                DO jk = 2 , imld(ji,jj)
727                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
728                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     &
729                        &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   &
730                        &                                          * zsc_uw_2(ji,jj)                                    )
731                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
732                END DO  ! jk loop
733             ELSE
734             ! stable conditions
735                DO jk = 2, ibld(ji,jj)
736                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
737                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
738                END DO
739             ENDIF
740          END DO        ! ji loop
741       END DO           ! jj loop
742
743       IF(ln_dia_osm) THEN
744          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
745          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
746       END IF
747! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
748
749       WHERE ( lconv )
750          zsc_wth_1 = zwth0
751          zsc_ws_1 = zws0
752       ELSEWHERE
753          zsc_wth_1 = 2.0 * zwthav
754          zsc_ws_1 = zws0
755       ENDWHERE
756
757       DO jj = 2, jpjm1
758          DO ji = 2, jpim1
759            IF ( lconv(ji,jj) ) THEN
760               DO jk = 2, imld(ji,jj)
761                  zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
762                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                &
763                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
764                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
765                       &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) )
766                  !
767                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  &
768                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
769                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
770                       &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) )
771               END DO
772            ELSE
773               DO jk = 2, ibld(ji,jj)
774                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
775                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
776                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
777               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
778                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
779               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
780               END DO
781            ENDIF
782          ENDDO    ! ji loop
783       END DO      ! jj loop
784
785       WHERE ( lconv )
786          zsc_uw_1 = zustar**2
787          zsc_vw_1 = ff_t * zustke * zhml
788       ELSEWHERE
789          zsc_uw_1 = zustar**2
790          zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
791          zsc_vw_1 = ff_t * zustke * zhbl
792          zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1
793       ENDWHERE
794
795       DO jj = 2, jpjm1
796          DO ji = 2, jpim1
797             IF ( lconv(ji,jj) ) THEN
798               DO jk = 2, imld(ji,jj)
799                  zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
800                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
801                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
802                       & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)
803                  !
804                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
805                       & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
806               END DO
807             ELSE
808               DO jk = 2, ibld(ji,jj)
809                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
810                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
811                  IF ( zznd_d <= 2.0 ) THEN
812                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
813                          &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
814                     !
815                  ELSE
816                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
817                          & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
818                     !
819                  ENDIF
820
821                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
822                       & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj)
823                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
824                       & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
825               END DO
826             ENDIF
827          END DO
828       END DO
829
830       IF(ln_dia_osm) THEN
831          IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
832          IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
833          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
834          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
835          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
836          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
837       END IF
838!
839! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
840
841      DO jj = 2, jpjm1
842         DO ji = 2, jpim1
843            IF ( lconv(ji,jj) ) THEN
844               DO jk = 2, ibld(ji,jj)
845                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
846                  IF ( znd >= 0.0 ) THEN
847                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
848                     ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
849                  ELSE
850                     ghamu(ji,jj,jk) = 0._wp
851                     ghamv(ji,jj,jk) = 0._wp
852                  ENDIF
853               END DO
854            ELSE
855               DO jk = 2, ibld(ji,jj)
856                  znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
857                  IF ( znd >= 0.0 ) THEN
858                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
859                     ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
860                  ELSE
861                     ghamu(ji,jj,jk) = 0._wp
862                     ghamv(ji,jj,jk) = 0._wp
863                  ENDIF
864               END DO
865            ENDIF
866         END DO
867      END DO
868
869       IF(ln_dia_osm) THEN
870          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
871          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
872       END IF
873      ! pynocline contributions
874       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small
875       zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln )
876       DO jj = 2, jpjm1
877          DO ji = 2, jpim1
878             IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
879                DO jk= 2, ibld(ji,jj)
880                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
881                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
882                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
883                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
884                   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)
885                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
886                END DO
887             END IF
888          END DO
889       END DO
890
891! Entrainment contribution.
892
893       DO jj=2, jpjm1
894          DO ji = 2, jpim1
895            IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN
896               DO jk = 1, imld(ji,jj) - 1
897                  znd=gdepw_n(ji,jj,jk) / zhml(ji,jj)
898                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
899                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
900                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd
901                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd
902               END DO
903               DO jk = imld(ji,jj), ibld(ji,jj)
904                  znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
905                  ! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
906                  ! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
907                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd )
908                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd )
909                END DO
910             ENDIF
911
912             ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
913             ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
914             ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
915             ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
916          END DO       ! ji loop
917       END DO          ! jj loop
918
919       IF(ln_dia_osm) THEN
920          IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
921          IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
922          IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse )
923          IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse )
924          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
925          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
926          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
927       END IF
928       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
929       ! Need to put in code for contributions that are applied explicitly to
930       ! the prognostic variables
931       !  1. Entrainment flux
932       !
933       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
934
935
936
937       ! rotate non-gradient velocity terms back to model reference frame
938
939       DO jj = 2, jpjm1
940          DO ji = 2, jpim1
941             DO jk = 2, ibld(ji,jj)
942                ztemp = ghamu(ji,jj,jk)
943                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
944                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
945             END DO
946          END DO
947       END DO
948
949       IF(ln_dia_osm) THEN
950          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
951          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
952          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
953       END IF
954
955! KPP-style Ri# mixing
956       IF( ln_kpprimix) THEN
957          DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
958             DO jj = 1, jpjm1
959                DO ji = 1, jpim1   ! vector opt.
960                   z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
961                        &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
962                        &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
963                   z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
964                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
965                        &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
966                END DO
967             END DO
968          END DO
969      !
970         DO jk = 2, jpkm1
971            DO jj = 2, jpjm1
972               DO ji = 2, jpim1   ! vector opt.
973                  !                                          ! shear prod. at w-point weightened by mask
974                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
975                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
976                  !                                          ! local Richardson number
977                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
978                  zfri =  MIN( zri / rn_riinfty , 1.0_wp )
979                  zfri  = ( 1.0_wp - zfri * zfri )
980                  zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk)
981                END DO
982             END DO
983          END DO
984
985          DO jj = 2, jpjm1
986             DO ji = 2, jpim1
987                DO jk = ibld(ji,jj) + 1, jpkm1
988                   zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
989                   zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
990                END DO
991             END DO
992          END DO
993
994       END IF ! ln_kpprimix = .true.
995
996! KPP-style set diffusivity large if unstable below BL
997       IF( ln_convmix) THEN
998          DO jj = 2, jpjm1
999             DO ji = 2, jpim1
1000                DO jk = ibld(ji,jj) + 1, jpkm1
1001                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
1002                END DO
1003             END DO
1004          END DO
1005       END IF ! ln_convmix = .true.
1006
1007       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
1008       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
1009
1010       ! GN 25/8: need to change tmask --> wmask
1011
1012     DO jk = 2, jpkm1
1013         DO jj = 2, jpjm1
1014             DO ji = 2, jpim1
1015                p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1016                p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
1017             END DO
1018         END DO
1019     END DO
1020      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid  (sign unchanged), needed to caclulate gham[uv] on u and v grids
1021     CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1.,   &
1022      &                  ghamu, 'W', 1. , ghamv, 'W', 1. )
1023       DO jk = 2, jpkm1
1024           DO jj = 2, jpjm1
1025               DO ji = 2, jpim1
1026                  ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
1027                     &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1028
1029                  ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
1030                      &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1031
1032                  ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk)
1033                  ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk)
1034               END DO
1035           END DO
1036        END DO
1037        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged)
1038        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged)
1039        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1.,   &
1040         &                  ghamu, 'U', -1. , ghamv, 'V', -1. )
1041
1042      IF(ln_dia_osm) THEN
1043         SELECT CASE (nn_osm_wave)
1044         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1).
1045         CASE(0:1)
1046            IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift
1047            IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift
1048            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
1049         ! Stokes drift read in from sbcwave  (=2).
1050         CASE(2)
1051            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift
1052            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift
1053            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                   ! wave mean period
1054            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                   ! significant wave height
1055            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
1056            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum
1057            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10
1058            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* &
1059                 & SQRT(ut0sd**2 + vt0sd**2 ) )
1060         END SELECT
1061         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )            ! <Tw_NL>
1062         IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams )            ! <Sw_NL>
1063         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )            ! <uw_NL>
1064         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )            ! <vw_NL>
1065         IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 )            ! <Tw_0>
1066         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0>
1067         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth
1068         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k
1069         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base
1070         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base
1071         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base
1072         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base
1073         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base
1074         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )               ! Initial boundary-layer depth
1075         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )               ! Initial boundary-layer depth
1076         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth
1077         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points
1078         IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc )         ! convective velocity scale
1079         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale
1080         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale
1081         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale
1082         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir #
1083         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
1084         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
1085         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine
1086         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine
1087         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine
1088         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine
1089         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine
1090         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux
1091         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux
1092         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux
1093         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux
1094         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML
1095      END IF
1096
1097CONTAINS
1098
1099
1100   ! Alan: do we need zb?
1101   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv )
1102     !!---------------------------------------------------------------------
1103     !!                   ***  ROUTINE zdf_vertical_average  ***
1104     !!
1105     !! ** Purpose : Determines vertical averages from surface to jnlev.
1106     !!
1107     !! ** Method  : Averages are calculated from the surface to jnlev.
1108     !!              The external level used to calculate differences is ibld+ibld_ext
1109     !!
1110     !!----------------------------------------------------------------------
1111
1112        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over.
1113
1114        ! Alan: do we need zb?
1115        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! Average temperature and salinity
1116        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components
1117        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
1118        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components.
1119
1120        INTEGER :: jk, ji, jj
1121        REAL(wp) :: zthick, zthermal, zbeta
1122
1123
1124        zt   = 0._wp
1125        zs   = 0._wp
1126        zu   = 0._wp
1127        zv   = 0._wp
1128        DO jj = 2, jpjm1                                 !  Vertical slab
1129         DO ji = 2, jpim1
1130            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1131            zbeta    = rab_n(ji,jj,1,jp_sal)
1132               ! average over depth of boundary layer
1133            zthick = epsln
1134            DO jk = 2, jnlev_av(ji,jj)
1135               zthick = zthick + e3t_n(ji,jj,jk)
1136               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
1137               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
1138               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) &
1139                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
1140                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
1141               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) &
1142                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
1143                     &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
1144            END DO
1145            zt(ji,jj) = zt(ji,jj) / zthick
1146            zs(ji,jj) = zs(ji,jj) / zthick
1147            zu(ji,jj) = zu(ji,jj) / zthick
1148            zv(ji,jj) = zv(ji,jj) / zthick
1149            ! Alan: do we need zb?
1150            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem)
1151            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal)
1152            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) &
1153                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) )
1154            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) &
1155                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) )
1156            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
1157         END DO
1158        END DO
1159   END SUBROUTINE zdf_osm_vertical_average
1160
1161   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
1162     !!---------------------------------------------------------------------
1163     !!                   ***  ROUTINE zdf_velocity_rotation  ***
1164     !!
1165     !! ** Purpose : Rotates frame of reference of averaged velocity components.
1166     !!
1167     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w
1168     !!
1169     !!----------------------------------------------------------------------
1170
1171        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle
1172        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current
1173        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline
1174
1175        INTEGER :: ji, jj
1176        REAL(wp) :: ztemp
1177
1178        DO jj = 2, jpjm1
1179           DO ji = 2, jpim1
1180              ztemp = zu(ji,jj)
1181              zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
1182              zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
1183              ztemp = zdu(ji,jj)
1184              zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
1185              zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj)
1186           END DO
1187        END DO
1188    END SUBROUTINE zdf_osm_velocity_rotation
1189
1190    SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz )
1191     !!---------------------------------------------------------------------
1192     !!                   ***  ROUTINE zdf_osm_external_gradients  ***
1193     !!
1194     !! ** Purpose : Calculates the gradients below the OSBL
1195     !!
1196     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient.
1197     !!
1198     !!----------------------------------------------------------------------
1199
1200     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy.
1201
1202     INTEGER :: jj, ji, jkb, jkb1
1203     REAL(wp) :: zthermal, zbeta
1204
1205
1206     DO jj = 2, jpjm1
1207        DO ji = 2, jpim1
1208           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1209              zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1210              zbeta    = rab_n(ji,jj,1,jp_sal)
1211              jkb = ibld(ji,jj) + ibld_ext
1212              jkb1 = MIN(jkb + 1, mbkt(ji,jj))
1213              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
1214                   &   / e3t_n(ji,jj,ibld(ji,jj))
1215              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
1216                   &   / e3t_n(ji,jj,ibld(ji,jj))
1217              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
1218           END IF
1219        END DO
1220     END DO
1221    END SUBROUTINE zdf_osm_external_gradients
1222
1223    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz )
1224
1225     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline
1226
1227     INTEGER :: jk, jj, ji
1228     REAL(wp) :: ztgrad, zsgrad, zbgrad
1229     REAL(wp) :: zgamma_b_nd, zgamma_c, znd
1230     REAL(wp) :: zzeta_s=0.3
1231
1232     DO jj = 2, jpjm1
1233        DO ji = 2, jpim1
1234           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1235              IF ( lconv(ji,jj) ) THEN  ! convective conditions
1236                 IF ( zdb_bl(ji,jj) > 0._wp .AND. zdbdz_ext(ji,jj) > 0._wp ) THEN
1237                    ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj)
1238                    zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj)
1239                    zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj)
1240                    zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1241                    zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
1242                         & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
1243                    DO jk = 2, ibld(ji,jj)+ibld_ext
1244                       znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1245                       IF ( znd <= zzeta_s ) THEN
1246                          zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * &
1247                               & EXP( -6.0 * ( znd -zzeta_s )**2 )
1248                          zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * &
1249                               & EXP( -6.0 * ( znd -zzeta_s )**2 )
1250                          zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * &
1251                               & EXP( -6.0 * ( znd -zzeta_s )**2 )
1252                       ELSE
1253                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1254                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1255                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1256                       ENDIF
1257                    END DO
1258                 ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing
1259              ELSE
1260                 ! stable conditions
1261                 ! if pycnocline profile only defined when depth steady of increasing.
1262                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady.
1263                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1264                       IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
1265                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
1266                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
1267                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
1268                          DO jk = 2, ibld(ji,jj)
1269                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1270                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1271                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1272                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1273                          END DO
1274                       ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
1275                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
1276                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
1277                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
1278                          DO jk = 2, ibld(ji,jj)
1279                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
1280                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1281                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1282                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1283                          END DO
1284                       ENDIF ! IF (zhol >=0.5)
1285                    ENDIF    ! IF (zdb_bl> 0.)
1286                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
1287              ENDIF          ! IF (lconv)
1288           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1289        END DO
1290     END DO
1291
1292    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
1293
1294    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
1295      !!---------------------------------------------------------------------
1296      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  ***
1297      !!
1298      !! ** Purpose : Calculates velocity shear in the pycnocline
1299      !!
1300      !! ** Method  :
1301      !!
1302      !!----------------------------------------------------------------------
1303
1304      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
1305
1306      INTEGER :: jk, jj, ji
1307      REAL(wp) :: zugrad, zvgrad, znd
1308      REAL(wp) :: zzeta_v = 0.45
1309      !
1310      DO jj = 2, jpjm1
1311         DO ji = 2, jpim1
1312            !
1313            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1314               IF ( lconv (ji,jj) ) THEN
1315                  ! Unstable conditions
1316                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
1317                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
1318                       &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
1319                  !Alan is this right?
1320                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
1321                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
1322                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
1323                       &      )/ (zdh(ji,jj)  + epsln )
1324                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
1325                     znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
1326                     IF ( znd <= 0.0 ) THEN
1327                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
1328                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
1329                     ELSE
1330                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
1331                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
1332                     ENDIF
1333                  END DO
1334               ELSE
1335                  ! stable conditions
1336                  zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
1337                  zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
1338                  DO jk = 2, ibld(ji,jj)
1339                     znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1340                     IF ( znd < 1.0 ) THEN
1341                        zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
1342                     ELSE
1343                        zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
1344                     ENDIF
1345                     zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
1346                  END DO
1347               ENDIF
1348               !
1349            END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1350         END DO
1351      END DO
1352    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
1353
1354    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
1355     !!---------------------------------------------------------------------
1356     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
1357     !!
1358     !! ** Purpose : Calculates the rate at which hbl changes.
1359     !!
1360     !! ** Method  :
1361     !!
1362     !!----------------------------------------------------------------------
1363
1364    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl
1365
1366    INTEGER :: jj, ji
1367    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert
1368    REAL(wp) :: zvel_max, zwb_min
1369    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
1370    REAL(wp) :: zzeta_m = 0.3
1371    REAL(wp) :: zgamma_c = 2.0
1372    REAL(wp) :: zdhoh = 0.1
1373    REAL(wp) :: alpha_bc = 0.5
1374
1375    DO jj = 2, jpjm1
1376       DO ji = 2, jpim1
1377          IF ( lconv(ji,jj) ) THEN    ! Convective
1378             ! Alan is this right?  Yes, it's a bit different from the previous relationship
1379             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) &
1380             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1381             zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
1382             zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
1383             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
1384             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
1385             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
1386                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
1387
1388             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
1389                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
1390                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
1391                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1392             !
1393             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj)
1394                  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) / &
1395                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1396                  zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1397                ! added ajgn 23 July as temporay fix
1398             zdhdt_2(ji,jj) = 0._wp
1399
1400                ! commented out ajgn 23 July as temporay fix
1401!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1402! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary.
1403!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1404!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj)
1405!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1406!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj)
1407!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min )
1408!                     ! Alan no idea what this should be?
1409!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) &
1410!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) &
1411!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj)
1412!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
1413!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN
1414!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj)
1415!                 ENDIF
1416            ELSE                        ! Stable
1417                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
1418                zdhdt_2(ji,jj) = 0._wp
1419                IF ( zdhdt(ji,jj) < 0._wp ) THEN
1420                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
1421                    zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
1422                ELSE
1423                    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) )
1424                ENDIF
1425                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert
1426            ENDIF
1427         END DO
1428      END DO
1429    END SUBROUTINE zdf_osm_calculate_dhdt
1430
1431    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
1432     !!---------------------------------------------------------------------
1433     !!                   ***  ROUTINE zdf_osm_timestep_hbl  ***
1434     !!
1435     !! ** Purpose : Increments hbl.
1436     !!
1437     !! ** Method  : If thechange in hbl exceeds one model level the change is
1438     !!              is calculated by moving down the grid, changing the buoyancy
1439     !!              jump. This is to ensure that the change in hbl does not
1440     !!              overshoot a stable layer.
1441     !!
1442     !!----------------------------------------------------------------------
1443
1444
1445    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2   ! rates of change of hbl.
1446
1447    INTEGER :: jk, jj, ji, jm
1448    REAL(wp) :: zhbl_s, zvel_max, zdb
1449    REAL(wp) :: zthermal, zbeta
1450
1451     DO jj = 2, jpjm1
1452         DO ji = 2, jpim1
1453           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
1454!
1455! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
1456!
1457              zhbl_s = hbl(ji,jj)
1458              jm = imld(ji,jj)
1459              zthermal = rab_n(ji,jj,1,jp_tem)
1460              zbeta = rab_n(ji,jj,1,jp_sal)
1461
1462
1463              IF ( lconv(ji,jj) ) THEN
1464!unstable
1465                    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) / &
1466                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1467                 DO jk = imld(ji,jj), ibld(ji,jj)
1468                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
1469                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
1470                         &  0.0 ) + zvel_max
1471
1472                      zhbl_s = zhbl_s + MIN( &
1473                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
1474                        & e3w_n(ji,jj,jm) )
1475                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1476
1477                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
1478                 END DO
1479                 hbl(ji,jj) = zhbl_s
1480                 ibld(ji,jj) = jm
1481             ELSE
1482! stable
1483                 DO jk = imld(ji,jj), ibld(ji,jj)
1484                    zdb = MAX( &
1485                         & grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&
1486                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
1487                         & 0.0 ) + &
1488             &       2.0 * zvstr(ji,jj)**2 / zhbl_s
1489
1490                    ! Alan is thuis right? I have simply changed hbli to hbl
1491                    zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
1492                    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) ) ) * &
1493               &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
1494                    zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
1495                    zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
1496
1497                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
1498                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
1499                 END DO
1500             ENDIF   ! IF ( lconv )
1501             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
1502             ibld(ji,jj) = MAX(jm, 4 )
1503           ELSE
1504! change zero or one model level.
1505             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
1506           ENDIF
1507           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
1508         END DO
1509      END DO
1510
1511    END SUBROUTINE zdf_osm_timestep_hbl
1512
1513    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
1514      !!---------------------------------------------------------------------
1515      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  ***
1516      !!
1517      !! ** Purpose : Calculates thickness of the pycnocline
1518      !!
1519      !! ** Method  : The thickness is calculated from a prognostic equation
1520      !!              that relaxes the pycnocine thickness to a diagnostic
1521      !!              value. The time change is calculated assuming the
1522      !!              thickness relaxes exponentially. This is done to deal
1523      !!              with large timesteps.
1524      !!
1525      !!----------------------------------------------------------------------
1526
1527      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness.
1528      !
1529      INTEGER :: jj, ji
1530      INTEGER :: inhml
1531      REAL(wp) :: zari, ztau, zddhdt
1532
1533
1534      DO jj = 2, jpjm1
1535         DO ji = 2, jpim1
1536            IF ( lconv(ji,jj) ) THEN
1537               IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1538                  IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
1539                     zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1540                          (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1541                  ELSE                                                     ! unstable
1542                     zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1543                          (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
1544                  ENDIF
1545                  ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
1546                  zddhdt = zari * hbl(ji,jj)
1547               ELSE
1548                  ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
1549                  zddhdt = 0.2 * hbl(ji,jj)
1550               ENDIF
1551               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1552               ! Alan: this hml is never defined or used
1553            ELSE   ! IF (lconv)
1554               ztau = hbl(ji,jj) / zvstr(ji,jj)
1555               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
1556                  ! boundary layer deepening
1557                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1558                     ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
1559                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
1560                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
1561                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj)
1562                  ELSE
1563                     zddhdt = 0.2 * hbl(ji,jj)
1564                  ENDIF
1565               ELSE     ! IF(dhdt < 0)
1566                  zddhdt = 0.2 * hbl(ji,jj)
1567               ENDIF    ! IF (dhdt >= 0)
1568               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
1569               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
1570               ! Alan: this hml is never defined or used -- do we need it?
1571            ENDIF       ! IF (lconv)
1572
1573            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
1574            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
1575            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
1576            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
1577            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1578         END DO
1579      END DO
1580
1581    END SUBROUTINE zdf_osm_pycnocline_thickness
1582
1583END SUBROUTINE zdf_osm
1584
1585
1586   SUBROUTINE zdf_osm_init
1587     !!----------------------------------------------------------------------
1588     !!                  ***  ROUTINE zdf_osm_init  ***
1589     !!
1590     !! ** Purpose :   Initialization of the vertical eddy diffivity and
1591     !!      viscosity when using a osm turbulent closure scheme
1592     !!
1593     !! ** Method  :   Read the namosm namelist and check the parameters
1594     !!      called at the first timestep (nit000)
1595     !!
1596     !! ** input   :   Namlist namosm
1597     !!----------------------------------------------------------------------
1598     INTEGER  ::   ios            ! local integer
1599     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1600     REAL z1_t2
1601     !!
1602     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
1603          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
1604          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave
1605     !!----------------------------------------------------------------------
1606     !
1607     REWIND( numnam_ref )              ! Namelist namzdf_osm in reference namelist : Osmosis ML model
1608     READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
1609901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist', lwp )
1610
1611     REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
1612     READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
1613902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist', lwp )
1614     IF(lwm) WRITE ( numond, namzdf_osm )
1615
1616     IF(lwp) THEN                    ! Control print
1617        WRITE(numout,*)
1618        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
1619        WRITE(numout,*) '~~~~~~~~~~~~'
1620        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters'
1621        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la
1622        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
1623        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
1624        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes
1625        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
1626        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave
1627        SELECT CASE (nn_osm_wave)
1628        CASE(0)
1629           WRITE(numout,*) '     calculated assuming constant La#=0.3'
1630        CASE(1)
1631           WRITE(numout,*) '     calculated from Pierson Moskowitz wind-waves'
1632        CASE(2)
1633           WRITE(numout,*) '     calculated from ECMWF wave fields'
1634        END SELECT
1635        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
1636        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix
1637        WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
1638        WRITE(numout,*) '     maximum shear diffusivity at Rig = 0    (m2/s) rn_difri = ', rn_difri
1639        WRITE(numout,*) '     Use large mixing below BL when unstable       ln_convmix = ', ln_convmix
1640        WRITE(numout,*) '     diffusivity when unstable below BL     (m2/s) rn_difconv = ', rn_difconv
1641     ENDIF
1642
1643     !                              ! allocate zdfosm arrays
1644     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
1645
1646
1647     call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle
1648
1649
1650     IF( ln_zdfddm) THEN
1651        IF(lwp) THEN
1652           WRITE(numout,*)
1653           WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
1654           WRITE(numout,*) '    CAUTION : done in routine zdfosm, not in routine zdfddm '
1655        ENDIF
1656     ENDIF
1657
1658
1659     !set constants not in namelist
1660     !-----------------------------
1661
1662     IF(lwp) THEN
1663        WRITE(numout,*)
1664     ENDIF
1665
1666     IF (nn_osm_wave == 0) THEN
1667        dstokes(:,:) = rn_osm_dstokes
1668     END IF
1669
1670     ! Horizontal average : initialization of weighting arrays
1671     ! -------------------
1672
1673     SELECT CASE ( nn_ave )
1674
1675     CASE ( 0 )                ! no horizontal average
1676        IF(lwp) WRITE(numout,*) '          no horizontal average on avt'
1677        IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
1678        ! weighting mean arrays etmean
1679        !           ( 1  1 )
1680        ! avt = 1/4 ( 1  1 )
1681        !
1682        etmean(:,:,:) = 0.e0
1683
1684        DO jk = 1, jpkm1
1685           DO jj = 2, jpjm1
1686              DO ji = 2, jpim1   ! vector opt.
1687                 etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
1688                      &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
1689                      &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
1690              END DO
1691           END DO
1692        END DO
1693
1694     CASE ( 1 )                ! horizontal average
1695        IF(lwp) WRITE(numout,*) '          horizontal average on avt'
1696        ! weighting mean arrays etmean
1697        !           ( 1/2  1  1/2 )
1698        ! avt = 1/8 ( 1    2  1   )
1699        !           ( 1/2  1  1/2 )
1700        etmean(:,:,:) = 0.e0
1701
1702        DO jk = 1, jpkm1
1703           DO jj = 2, jpjm1
1704              DO ji = 2, jpim1   ! vector opt.
1705                 etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
1706                      & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
1707                      &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
1708                      &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
1709                      &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
1710                      &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
1711              END DO
1712           END DO
1713        END DO
1714
1715     CASE DEFAULT
1716        WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
1717        CALL ctl_stop( ctmp1 )
1718
1719     END SELECT
1720
1721     ! Initialization of vertical eddy coef. to the background value
1722     ! -------------------------------------------------------------
1723     DO jk = 1, jpk
1724        avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1725     END DO
1726
1727     ! zero the surface flux for non local term and osm mixed layer depth
1728     ! ------------------------------------------------------------------
1729     ghamt(:,:,:) = 0.
1730     ghams(:,:,:) = 0.
1731     ghamu(:,:,:) = 0.
1732     ghamv(:,:,:) = 0.
1733     !
1734     IF( lwxios ) THEN
1735        CALL iom_set_rstw_var_active('wn')
1736        CALL iom_set_rstw_var_active('hbl')
1737        CALL iom_set_rstw_var_active('hbli')
1738     ENDIF
1739   END SUBROUTINE zdf_osm_init
1740
1741
1742   SUBROUTINE osm_rst( kt, cdrw )
1743     !!---------------------------------------------------------------------
1744     !!                   ***  ROUTINE osm_rst  ***
1745     !!
1746     !! ** Purpose :   Read or write BL fields in restart file
1747     !!
1748     !! ** Method  :   use of IOM library. If the restart does not contain
1749     !!                required fields, they are recomputed from stratification
1750     !!----------------------------------------------------------------------
1751
1752     INTEGER, INTENT(in) :: kt
1753     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1754
1755     INTEGER ::   id1, id2   ! iom enquiry index
1756     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1757     INTEGER  ::   iiki, ikt ! local integer
1758     REAL(wp) ::   zhbf           ! tempory scalars
1759     REAL(wp) ::   zN2_c           ! local scalar
1760     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
1761     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
1762     !!----------------------------------------------------------------------
1763     !
1764     !!-----------------------------------------------------------------------------
1765     ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
1766     !!-----------------------------------------------------------------------------
1767     IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
1768        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. )
1769        IF( id1 > 0 ) THEN                       ! 'wn' exists; read
1770           CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios )
1771           WRITE(numout,*) ' ===>>>> :  wn read from restart file'
1772        ELSE
1773           wn(:,:,:) = 0._wp
1774           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
1775        END IF
1776
1777        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
1778        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. )
1779        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return
1780           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
1781           CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios  )
1782           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file'
1783           RETURN
1784        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate
1785           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
1786        END IF
1787     END IF
1788
1789     !!-----------------------------------------------------------------------------
1790     ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
1791     !!-----------------------------------------------------------------------------
1792     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
1793        IF(lwp) WRITE(numout,*) '---- osm-rst ----'
1794         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios )
1795         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios )
1796         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh, ldxios = lwxios )
1797        RETURN
1798     END IF
1799
1800     !!-----------------------------------------------------------------------------
1801     ! Getting hbl, no restart file with hbl, so calculate from surface stratification
1802     !!-----------------------------------------------------------------------------
1803     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
1804     ! w-level of the mixing and mixed layers
1805     CALL eos_rab( tsn, rab_n )
1806     CALL bn2(tsn, rab_n, rn2)
1807     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
1808     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
1809     zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
1810     !
1811     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
1812     DO jk = 1, jpkm1
1813        DO jj = 1, jpj              ! Mixed layer level: w-level
1814           DO ji = 1, jpi
1815              ikt = mbkt(ji,jj)
1816              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
1817              IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
1818           END DO
1819        END DO
1820     END DO
1821     !
1822     DO jj = 1, jpj
1823        DO ji = 1, jpi
1824           iiki = MAX(4,imld_rst(ji,jj))
1825           hbl (ji,jj) = gdepw_n(ji,jj,iiki  )    ! Turbocline depth
1826           dh (ji,jj) = e3t_n(ji,jj,iiki-1  )     ! Turbocline depth
1827        END DO
1828     END DO
1829     WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
1830     wn(:,:,:) = 0._wp
1831     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
1832   END SUBROUTINE osm_rst
1833
1834
1835   SUBROUTINE tra_osm( kt )
1836      !!----------------------------------------------------------------------
1837      !!                  ***  ROUTINE tra_osm  ***
1838      !!
1839      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
1840      !!
1841      !! ** Method  :   ???
1842      !!----------------------------------------------------------------------
1843      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
1844      !!----------------------------------------------------------------------
1845      INTEGER, INTENT(in) :: kt
1846      INTEGER :: ji, jj, jk
1847      !
1848      IF( kt == nit000 ) THEN
1849         IF(lwp) WRITE(numout,*)
1850         IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
1851         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1852      ENDIF
1853
1854      IF( l_trdtra )   THEN                    !* Save ta and sa trends
1855         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
1856         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
1857      ENDIF
1858
1859      DO jk = 1, jpkm1
1860         DO jj = 2, jpjm1
1861            DO ji = 2, jpim1
1862               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
1863                  &                 - (  ghamt(ji,jj,jk  )  &
1864                  &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
1865               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
1866                  &                 - (  ghams(ji,jj,jk  )  &
1867                  &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
1868            END DO
1869         END DO
1870      END DO
1871
1872      ! save the non-local tracer flux trends for diagnostics
1873      IF( l_trdtra )   THEN
1874         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
1875         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
1876
1877         CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt )
1878         CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds )
1879         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
1880      ENDIF
1881
1882      IF(ln_ctl) THEN
1883         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   &
1884         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
1885      ENDIF
1886      !
1887   END SUBROUTINE tra_osm
1888
1889
1890   SUBROUTINE trc_osm( kt )          ! Dummy routine
1891      !!----------------------------------------------------------------------
1892      !!                  ***  ROUTINE trc_osm  ***
1893      !!
1894      !! ** Purpose :   compute and add to the passive tracer trend the non-local
1895      !!                 passive tracer flux
1896      !!
1897      !!
1898      !! ** Method  :   ???
1899      !!----------------------------------------------------------------------
1900      !
1901      !!----------------------------------------------------------------------
1902      INTEGER, INTENT(in) :: kt
1903      WRITE(*,*) 'trc_osm: Not written yet', kt
1904   END SUBROUTINE trc_osm
1905
1906
1907   SUBROUTINE dyn_osm( kt )
1908      !!----------------------------------------------------------------------
1909      !!                  ***  ROUTINE dyn_osm  ***
1910      !!
1911      !! ** Purpose :   compute and add to the velocity trend the non-local flux
1912      !! copied/modified from tra_osm
1913      !!
1914      !! ** Method  :   ???
1915      !!----------------------------------------------------------------------
1916      INTEGER, INTENT(in) ::   kt   !
1917      !
1918      INTEGER :: ji, jj, jk   ! dummy loop indices
1919      !!----------------------------------------------------------------------
1920      !
1921      IF( kt == nit000 ) THEN
1922         IF(lwp) WRITE(numout,*)
1923         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
1924         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1925      ENDIF
1926      !code saving tracer trends removed, replace with trdmxl_oce
1927
1928      DO jk = 1, jpkm1           ! add non-local u and v fluxes
1929         DO jj = 2, jpjm1
1930            DO ji = 2, jpim1
1931               ua(ji,jj,jk) =  ua(ji,jj,jk)                      &
1932                  &                 - (  ghamu(ji,jj,jk  )  &
1933                  &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
1934               va(ji,jj,jk) =  va(ji,jj,jk)                      &
1935                  &                 - (  ghamv(ji,jj,jk  )  &
1936                  &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
1937            END DO
1938         END DO
1939      END DO
1940      !
1941      ! code for saving tracer trends removed
1942      !
1943   END SUBROUTINE dyn_osm
1944
1945   !!======================================================================
1946
1947END MODULE zdfosm
Note: See TracBrowser for help on using the repository browser.