source: NEMO/trunk/src/OCE/ZDF/zdfosm.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 8 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

  • Property svn:keywords set to Id
File size: 86.7 KB
Line 
1MODULE zdfosm
2   !!======================================================================
3   !!                       ***  MODULE  zdfosm  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the OSMOSIS
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !!  History : NEMO 4.0  ! A. Grant, G. Nurser
8   !! 15/03/2017  Changed calculation of pycnocline thickness in unstable conditions and stable conditions AG
9   !! 15/03/2017  Calculation of pycnocline gradients for stable conditions changed. Pycnocline gradients now depend on stability of the OSBL. A.G
10   !! 06/06/2017  (1) Checks on sign of buoyancy jump in calculation of  OSBL depth.  A.G.
11   !!             (2) Removed variable zbrad0, zbradh and zbradav since they are not used.
12   !!             (3) Approximate treatment for shear turbulence.
13   !!                        Minimum values for zustar and zustke.
14   !!                        Add velocity scale, zvstr, that tends to zustar for large Langmuir numbers.
15   !!                        Limit maximum value for Langmuir number.
16   !!                        Use zvstr in definition of stability parameter zhol.
17   !!             (4) Modified parametrization of entrainment flux, changing original coefficient 0.0485 for Langmuir contribution to 0.135 * zla
18   !!             (5) For stable boundary layer add factor that depends on length of timestep to 'slow' collapse and growth. Make sure buoyancy jump not negative.
19   !!             (6) For unstable conditions when growth is over multiple levels, limit change to maximum of one level per cycle through loop.
20   !!             (7) Change lower limits for loops that calculate OSBL averages from 1 to 2. Large gradients between levels 1 and 2 can cause problems.
21   !!             (8) Change upper limits from ibld-1 to ibld.
22   !!             (9) Calculation of pycnocline thickness in unstable conditions. Check added to ensure that buoyancy jump is positive before calculating Ri.
23   !!            (10) Thickness of interface layer at base of the stable OSBL set by Richardson number. Gives continuity in transition from unstable OSBL.
24   !!            (11) Checks that buoyancy jump is poitive when calculating pycnocline profiles.
25   !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity.
26   !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information
27   !!            (14) Bouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out).
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   !!----------------------------------------------------------------------
32
33   !!----------------------------------------------------------------------
34   !!   'ln_zdfosm'                                             OSMOSIS scheme
35   !!----------------------------------------------------------------------
36   !!   zdf_osm       : update momentum and tracer Kz from osm scheme
37   !!   zdf_osm_init  : initialization, namelist read, and parameters control
38   !!   osm_rst       : read (or initialize) and write osmosis restart fields
39   !!   tra_osm       : compute and add to the T & S trend the non-local flux
40   !!   trc_osm       : compute and add to the passive tracer trend the non-local flux (TBD)
41   !!   dyn_osm       : compute and add to u & v trensd the non-local flux
42   !!----------------------------------------------------------------------
43   USE oce            ! ocean dynamics and active tracers
44                      ! uses ww from previous time step (which is now wb) to calculate hbl
45   USE dom_oce        ! ocean space and time domain
46   USE zdf_oce        ! ocean vertical physics
47   USE sbc_oce        ! surface boundary condition: ocean
48   USE sbcwave        ! surface wave parameters
49   USE phycst         ! physical constants
50   USE eosbn2         ! equation of state
51   USE traqsr         ! details of solar radiation absorption
52   USE zdfddm         ! double diffusion mixing (avs array)
53   USE iom            ! I/O library
54   USE lib_mpp        ! MPP library
55   USE trd_oce        ! ocean trends definition
56   USE trdtra         ! tracers trends
57   !
58   USE in_out_manager ! I/O manager
59   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
60   USE prtctl         ! Print control
61   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
62
63   IMPLICIT NONE
64   PRIVATE
65
66   PUBLIC   zdf_osm       ! routine called by step.F90
67   PUBLIC   zdf_osm_init  ! routine called by nemogcm.F90
68   PUBLIC   osm_rst       ! routine called by step.F90
69   PUBLIC   tra_osm       ! routine called by step.F90
70   PUBLIC   trc_osm       ! routine called by trcstp.F90
71   PUBLIC   dyn_osm       ! routine called by 'step.F90'
72
73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux
74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv    !: non-local v-momentum flux
75   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamt    !: non-local temperature flux (gamma/<ws>o)
76   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghams    !: non-local salinity flux (gamma/<ws>o)
77   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt
78   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth
79   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbli     !: intial boundary layer depth for stable blayer
80   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift.
81
82   !                      !!** Namelist  namzdf_osm  **
83   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la
84   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number
85   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift
86   REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs
87   INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt
88   INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
89   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la
90
91
92   LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing
93   REAL(wp) ::   rn_riinfty   = 0.7     ! local Richardson Number limit for shear instability
94   REAL(wp) ::   rn_difri    =  0.005   ! maximum shear mixing at Rig = 0    (m2/s)
95   LOGICAL  ::   ln_convmix  = .true.   ! Convective instability mixing
96   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s)
97
98   !                                    !!! ** General constants  **
99   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number
100   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3
101   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3
102
103   INTEGER :: idebug = 236
104   INTEGER :: jdebug = 228
105   !! * Substitutions
106#  include "do_loop_substitute.h90"
107   !!----------------------------------------------------------------------
108   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
109   !! $Id$
110   !! Software governed by the CeCILL license (see ./LICENSE)
111   !!----------------------------------------------------------------------
112CONTAINS
113
114   INTEGER FUNCTION zdf_osm_alloc()
115      !!----------------------------------------------------------------------
116      !!                 ***  FUNCTION zdf_osm_alloc  ***
117      !!----------------------------------------------------------------------
118     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), &
119          &      hbl(jpi,jpj)    ,  hbli(jpi,jpj)    , dstokes(jpi, jpj) ,                     &
120          &   etmean(jpi,jpj,jpk),  STAT= zdf_osm_alloc )
121     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
122     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
123   END FUNCTION zdf_osm_alloc
124
125
126   SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt )
127      !!----------------------------------------------------------------------
128      !!                   ***  ROUTINE zdf_osm  ***
129      !!
130      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
131      !!      coefficients and non local mixing using the OSMOSIS scheme
132      !!
133      !! ** Method :   The boundary layer depth hosm is diagnosed at tracer points
134      !!      from profiles of buoyancy, and shear, and the surface forcing.
135      !!      Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
136      !!
137      !!                      Kx =  hosm  Wx(sigma) G(sigma)
138      !!
139      !!             and the non local term ghamt = Cs / Ws(sigma) / hosm
140      !!      Below hosm  the coefficients are the sum of mixing due to internal waves
141      !!      shear instability and double diffusion.
142      !!
143      !!      -1- Compute the now interior vertical mixing coefficients at all depths.
144      !!      -2- Diagnose the boundary layer depth.
145      !!      -3- Compute the now boundary layer vertical mixing coefficients.
146      !!      -4- Compute the now vertical eddy vicosity and diffusivity.
147      !!      -5- Smoothing
148      !!
149      !!        N.B. The computation is done from jk=2 to jpkm1
150      !!             Surface value of avt are set once a time to zero
151      !!             in routine zdf_osm_init.
152      !!
153      !! ** Action  :   update the non-local terms ghamts
154      !!                update avt (before vertical eddy coef.)
155      !!
156      !! References : Large W.G., Mc Williams J.C. and Doney S.C.
157      !!         Reviews of Geophysics, 32, 4, November 1994
158      !!         Comments in the code refer to this paper, particularly
159      !!         the equation number. (LMD94, here after)
160      !!----------------------------------------------------------------------
161      INTEGER                   , INTENT(in   ) ::  kt             ! ocean time step
162      INTEGER                   , INTENT(in   ) ::  Kbb, Kmm, Krhs ! ocean time level indices
163      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points)
164      !!
165      INTEGER ::   ji, jj, jk                   ! dummy loop indices
166      INTEGER ::   ikbot, jkmax, jkm1, jkp2     !
167
168      REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube      !
169      REAL(wp) ::   zbeta, zthermal                                  !
170      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
171      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm       !
172      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density
173      INTEGER  ::   jm                          ! dummy loop indices
174      REAL(wp) ::   zr1, zr2, zr3, zr4, zrhop   ! Compression terms
175      REAL(wp) ::   zflag, zrn2, zdep21, zdep32, zdep43
176      REAL(wp) ::   zesh2, zri, zfri            ! Interior richardson mixing
177      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
178      REAL(wp) :: zt,zs,zu,zv,zrh               ! variables used in constructing averages
179! Scales
180      REAL(wp), DIMENSION(jpi,jpj) :: zrad0     ! Surface solar temperature flux (deg m/s)
181      REAL(wp), DIMENSION(jpi,jpj) :: zradh     ! Radiative flux at bl base (Buoyancy units)
182      REAL(wp), DIMENSION(jpi,jpj) :: zradav    ! Radiative flux, bl average (Buoyancy Units)
183      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
184      REAL(wp), DIMENSION(jpi,jpj) :: zwstrl    ! Langmuir velocity scale
185      REAL(wp), DIMENSION(jpi,jpj) :: zvstr     ! Velocity scale that ends to zustar for large Langmuir number.
186      REAL(wp), DIMENSION(jpi,jpj) :: zwstrc    ! Convective velocity scale
187      REAL(wp), DIMENSION(jpi,jpj) :: zuw0      ! Surface u-momentum flux
188      REAL(wp), DIMENSION(jpi,jpj) :: zvw0      ! Surface v-momentum flux
189      REAL(wp), DIMENSION(jpi,jpj) :: zwth0     ! Surface heat flux (Kinematic)
190      REAL(wp), DIMENSION(jpi,jpj) :: zws0      ! Surface freshwater flux
191      REAL(wp), DIMENSION(jpi,jpj) :: zwb0      ! Surface buoyancy flux
192      REAL(wp), DIMENSION(jpi,jpj) :: zwthav    ! Heat flux - bl average
193      REAL(wp), DIMENSION(jpi,jpj) :: zwsav     ! freshwater flux - bl average
194      REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average
195      REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux
196      REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift
197      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number
198      REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
199      REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
200      REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer
201      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl
202
203      ! mixed-layer variables
204
205      INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
206      INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
207
208      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
209      REAL(wp) :: zugrad,zvgrad        ! temporary variables for calculating pycnocline shear
210
211      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid
212      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid
213      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid
214      REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
215      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer
216      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer
217      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
218      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
219      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
220      REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline
221      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline
222      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline
223      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline
224      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc    ! u-shear across the pycnocline
225      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc    ! v-shear across the pycnocline
226
227      ! Flux-gradient relationship variables
228
229      REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale.
230
231      REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity
232      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
233      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.
234      REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
235
236      ! For calculating Ri#-dependent mixing
237      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du   ! u-shear^2
238      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv   ! v-shear^2
239      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
240
241      ! Temporary variables
242      INTEGER :: inhml
243      INTEGER :: i_lconv_alloc
244      REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
245      REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables
246      REAL(wp) :: zthick, zz0, zz1 ! temporary variables
247      REAL(wp) :: zvel_max, zhbl_s ! temporary variables
248      REAL(wp) :: zfac             ! temporary variable
249      REAL(wp) :: zus_x, zus_y     ! temporary Stokes drift
250      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
251      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
252
253      ! For debugging
254      INTEGER :: ikt
255      !!--------------------------------------------------------------------
256      !
257      ALLOCATE( lconv(jpi,jpj),  STAT= i_lconv_alloc )
258      IF( i_lconv_alloc /= 0 )   CALL ctl_warn('zdf_osm: failed to allocate lconv')
259
260      ibld(:,:)   = 0     ; imld(:,:)  = 0
261      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp
262      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp
263      zvw0(:,:)   = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:)      = 0._wp ; zwb0(:,:)      = 0._wp
264      zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:)     = 0._wp ; zwb_ent(:,:)   = 0._wp
265      zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
266      zhol(:,:)   = 0._wp
267      lconv(:,:)  = .FALSE.
268      ! mixed layer
269      ! no initialization of zhbl or zhml (or zdh?)
270      zhbl(:,:)    = 1._wp ; zhml(:,:)    = 1._wp ; zdh(:,:)      = 1._wp ; zdhdt(:,:)   = 0._wp
271      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp
272      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp
273      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp
274      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp
275      zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp
276      zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:)  = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp
277      zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp
278      !
279      zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
280      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
281      !
282      ! Flux-Gradient arrays.
283      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp
284      zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp
285      zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp
286      zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp
287      zhbl_t(:,:)     = 0._wp ; zdhdt(:,:)      = 0._wp
288
289      zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
290      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp
291
292      ! hbl = MAX(hbl,epsln)
293      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
294      ! Calculate boundary layer scales
295      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
296
297      ! Assume two-band radiation model for depth of OSBL
298     zz0 =       rn_abs       ! surface equi-partition in 2-bands
299     zz1 =  1. - rn_abs
300     DO_2D_00_00
301        ! Surface downward irradiance (so always +ve)
302        zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp
303        ! Downwards irradiance at base of boundary layer
304        zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
305        ! Downwards irradiance averaged over depth of the OSBL
306        zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
307              &                         + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
308     END_2D
309     ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
310     DO_2D_00_00
311        zthermal = rab_n(ji,jj,1,jp_tem)
312        zbeta    = rab_n(ji,jj,1,jp_sal)
313        ! Upwards surface Temperature flux for non-local term
314        zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1)
315        ! Upwards surface salinity flux for non-local term
316        zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)
317        ! Non radiative upwards surface buoyancy flux
318        zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj)
319        ! turbulent heat flux averaged over depth of OSBL
320        zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
321        ! turbulent salinity flux averaged over depth of the OBSL
322        zwsav(ji,jj) = 0.5 * zws0(ji,jj)
323        ! turbulent buoyancy flux averaged over the depth of the OBSBL
324        zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj)
325        ! Surface upward velocity fluxes
326        zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
327        zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
328        ! Friction velocity (zustar), at T-point : LMD94 eq. 2
329        zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
330        zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
331        zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
332     END_2D
333     ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
334     SELECT CASE (nn_osm_wave)
335     ! Assume constant La#=0.3
336     CASE(0)
337        DO_2D_00_00
338           zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
339           zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
340           zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 )
341           ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
342        END_2D
343     ! Assume Pierson-Moskovitz wind-wave spectrum
344     CASE(1)
345        DO_2D_00_00
346           ! Use wind speed wndm included in sbc_oce module
347           zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
348           dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav
349        END_2D
350     ! Use ECMWF wave fields as output from SBCWAVE
351     CASE(2)
352        zfac =  2.0_wp * rpi / 16.0_wp
353        DO_2D_00_00
354           ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
355           !    The coefficient 0.8 gives La=0.3  in this situation.
356           ! It could represent the effects of the spread of wave directions
357           ! around the mean wind. The effect of this adjustment needs to be tested.
358           zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), &
359                &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  )
360           dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
361        END_2D
362     END SELECT
363
364     ! Langmuir velocity scale (zwstrl), La # (zla)
365     ! mixed scale (zvstr), convective velocity scale (zwstrc)
366     DO_2D_00_00
367        ! Langmuir velocity scale (zwstrl), at T-point
368        zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
369        ! Modify zwstrl to allow for small and large values of dstokes/hbl.
370        ! Intended as a possible test. Doesn't affect LES results for entrainment,
371        !  but hasn't been shown to be correct as dstokes/h becomes large or small.
372        zwstrl(ji,jj) = zwstrl(ji,jj) *  &
373             & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * &
374             & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) ))
375        ! define La this way so effects of Stokes penetration depth on velocity scale are included
376        zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3
377        ! Velocity scale that tends to zustar for large Langmuir numbers
378        zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + &
379             & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
380
381        ! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
382        ! Note zustke and zwstrl are not amended.
383        IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45
384        !
385        ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
386        IF ( zwbav(ji,jj) > 0.0) THEN
387           zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
388           zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
389           lconv(ji,jj) = .TRUE.
390        ELSE
391           zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln )
392           lconv(ji,jj) = .FALSE.
393        ENDIF
394     END_2D
395
396     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
397     ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
398     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
399     ! BL must be always 2 levels deep.
400      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) )
401      ibld(:,:) = 3
402      DO_3D_00_00( 4, jpkm1 )
403         IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
404            ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
405         ENDIF
406      END_3D
407
408      DO_2D_00_00
409            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
410            zbeta    = rab_n(ji,jj,1,jp_sal)
411            zt   = 0._wp
412            zs   = 0._wp
413            zu   = 0._wp
414            zv   = 0._wp
415            ! average over depth of boundary layer
416            zthick=0._wp
417            DO jm = 2, ibld(ji,jj)
418               zthick=zthick+e3t(ji,jj,jm,Kmm)
419               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
420               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
421               zu   = zu  + e3t(ji,jj,jm,Kmm) &
422                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
423                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
424               zv   = zv  + e3t(ji,jj,jm,Kmm) &
425                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
426                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
427            END DO
428            zt_bl(ji,jj) = zt / zthick
429            zs_bl(ji,jj) = zs / zthick
430            zu_bl(ji,jj) = zu / zthick
431            zv_bl(ji,jj) = zv / zthick
432            zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
433            zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
434            zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
435                  &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
436            zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
437                  &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
438            zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
439            IF ( lconv(ji,jj) ) THEN    ! Convective
440                   zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) &
441                        &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) )
442
443                   zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) &
444                        &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
445! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment.
446!                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) &
447!                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) )
448
449!                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / &
450!                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
451                   zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) )
452            ELSE                        ! Stable
453                   zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) &
454                        &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) &
455                        & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) &
456                        &  * zwstrl(ji,jj)**3 / hbli(ji,jj)
457                   zzdhdt = zzdhdt + zwbav(ji,jj)
458                   IF ( zzdhdt < 0._wp ) THEN
459                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
460                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj)
461                   ELSE
462                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) &
463                           &  + MAX( zdb_bl(ji,jj), 0.0 )
464                   ENDIF
465                   zzdhdt = 2.0 * zzdhdt / zpert
466            ENDIF
467            zdhdt(ji,jj) = zzdhdt
468      END_2D
469
470      ! Calculate averages over depth of boundary layer
471      imld = ibld           ! use imld to hold previous blayer index
472      ibld(:,:) = 3
473
474      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it
475      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:))
476      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom
477
478      DO_3D_00_00( 4, jpkm1 )
479         IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
480            ibld(ji,jj) =  MIN(mbkt(ji,jj), jk)
481         ENDIF
482      END_3D
483
484!
485! Step through model levels taking account of buoyancy change to determine the effect on dhdt
486!
487      DO_2D_00_00
488         IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
489!
490! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
491!
492            zhbl_s = hbl(ji,jj)
493            jm = imld(ji,jj)
494            zthermal = rab_n(ji,jj,1,jp_tem)
495            zbeta = rab_n(ji,jj,1,jp_sal)
496            IF ( lconv(ji,jj) ) THEN
497!unstable
498               zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) &
499                    &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
500
501               DO jk = imld(ji,jj), ibld(ji,jj)
502                  zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) &
503                       & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max
504
505                  zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) )
506                  zhbl_s = MIN(zhbl_s, ht(ji,jj))
507
508                  IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1
509               END DO
510               hbl(ji,jj) = zhbl_s
511               ibld(ji,jj) = jm
512               hbli(ji,jj) = hbl(ji,jj)
513            ELSE
514! stable
515               DO jk = imld(ji,jj), ibld(ji,jj)
516                  zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )          &
517                       &               - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) &
518                       & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s
519
520                  zhbl_s = zhbl_s +  (                                                                                &
521                       &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     &
522                       &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               &
523                       &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   &
524                       &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) &
525                       &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include ww here
526
527                  zhbl_s = MIN(zhbl_s, ht(ji,jj))
528                  IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1
529               END DO
530               hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) )
531               ibld(ji,jj) = MAX(jm, 3 )
532               IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
533            ENDIF   ! IF ( lconv )
534         ELSE
535! change zero or one model level.
536            hbl(ji,jj) = zhbl_t(ji,jj)
537            IF ( lconv(ji,jj) ) THEN
538               hbli(ji,jj) = hbl(ji,jj)
539            ELSE
540               hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) )
541               IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
542            ENDIF
543         ENDIF
544         zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
545      END_2D
546      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms.
547
548! Recalculate averages over boundary layer after depth updated
549     ! Consider later  combining this into the loop above and looking for columns
550     ! where the index for base of the boundary layer have changed
551      DO_2D_00_00
552            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
553            zbeta    = rab_n(ji,jj,1,jp_sal)
554            zt   = 0._wp
555            zs   = 0._wp
556            zu   = 0._wp
557            zv   = 0._wp
558            ! average over depth of boundary layer
559            zthick=0._wp
560            DO jm = 2, ibld(ji,jj)
561               zthick=zthick+e3t(ji,jj,jm,Kmm)
562               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
563               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
564               zu   = zu  + e3t(ji,jj,jm,Kmm) &
565                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
566                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
567               zv   = zv  + e3t(ji,jj,jm,Kmm) &
568                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
569                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
570            END DO
571            zt_bl(ji,jj) = zt / zthick
572            zs_bl(ji,jj) = zs / zthick
573            zu_bl(ji,jj) = zu / zthick
574            zv_bl(ji,jj) = zv / zthick
575            zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
576            zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
577            zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
578                   &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
579            zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
580                   &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
581            zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
582            zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
583            IF ( lconv(ji,jj) ) THEN
584               IF ( zdb_bl(ji,jj) > 0._wp )THEN
585                  IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
586                        zari = 4.5 * ( zvstr(ji,jj)**2 ) &
587                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01
588                  ELSE                                                     ! unstable
589                        zari = 4.5 * ( zwstrc(ji,jj)**2 ) &
590                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01
591                  ENDIF
592                  IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline
593                     zari = 0.2
594                     zwb_ent(ji,jj) = 0._wp
595                  ENDIF
596                  inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 )
597                  imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
598                  zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
599                  zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
600               ELSE  ! IF (zdb_bl)
601                  imld(ji,jj) = ibld(ji,jj) - 1
602                  zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
603                  zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
604               ENDIF
605            ELSE   ! IF (lconv)
606               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
607               ! boundary layer deepening
608                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN
609               ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
610                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
611                       & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
612                     inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 )
613                     imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
614                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
615                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
616                  ELSE
617                     imld(ji,jj) = ibld(ji,jj) - 1
618                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
619                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
620                  ENDIF ! IF (zdb_bl > 0.0)
621               ELSE     ! IF(dhdt >= 0)
622               ! boundary layer collapsing.
623                  imld(ji,jj) = ibld(ji,jj)
624                  zhml(ji,jj) = zhbl(ji,jj)
625                  zdh(ji,jj) = 0._wp
626               ENDIF    ! IF (dhdt >= 0)
627            ENDIF       ! IF (lconv)
628      END_2D
629
630      ! Average over the depth of the mixed layer in the convective boundary layer
631      ! Also calculate entrainment fluxes for temperature and salinity
632      DO_2D_00_00
633         zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
634         zbeta    = rab_n(ji,jj,1,jp_sal)
635         IF ( lconv(ji,jj) ) THEN
636            zt   = 0._wp
637            zs   = 0._wp
638            zu   = 0._wp
639            zv   = 0._wp
640            ! average over depth of boundary layer
641            zthick=0._wp
642            DO jm = 2, imld(ji,jj)
643               zthick=zthick+e3t(ji,jj,jm,Kmm)
644               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
645               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
646               zu   = zu  + e3t(ji,jj,jm,Kmm) &
647                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
648                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
649               zv   = zv  + e3t(ji,jj,jm,Kmm) &
650                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
651                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
652            END DO
653            zt_ml(ji,jj) = zt / zthick
654            zs_ml(ji,jj) = zs / zthick
655            zu_ml(ji,jj) = zu / zthick
656            zv_ml(ji,jj) = zv / zthick
657            zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
658            zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
659            zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
660                  &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
661            zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
662                  &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
663            zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
664         ELSE
665         ! stable, if entraining calulate average below interface layer.
666            IF ( zdhdt(ji,jj) >= 0._wp ) THEN
667               zt   = 0._wp
668               zs   = 0._wp
669               zu   = 0._wp
670               zv   = 0._wp
671               ! average over depth of boundary layer
672               zthick=0._wp
673               DO jm = 2, imld(ji,jj)
674                  zthick=zthick+e3t(ji,jj,jm,Kmm)
675                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
676                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
677                  zu   = zu  + e3t(ji,jj,jm,Kmm) &
678                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
679                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
680                  zv   = zv  + e3t(ji,jj,jm,Kmm) &
681                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
682                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
683               END DO
684               zt_ml(ji,jj) = zt / zthick
685               zs_ml(ji,jj) = zs / zthick
686               zu_ml(ji,jj) = zu / zthick
687               zv_ml(ji,jj) = zv / zthick
688               zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
689               zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
690               zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
691                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
692               zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
693                     &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
694               zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
695            ENDIF
696         ENDIF
697      END_2D
698    !
699    ! rotate mean currents and changes onto wind align co-ordinates
700    !
701
702      DO_2D_00_00
703         ztemp = zu_ml(ji,jj)
704         zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj)
705         zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
706         ztemp = zdu_ml(ji,jj)
707         zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj)
708         zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
709 !
710         ztemp = zu_bl(ji,jj)
711         zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj)
712         zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
713         ztemp = zdu_bl(ji,jj)
714         zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj)
715         zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
716      END_2D
717
718     zuw_bse = 0._wp
719     zvw_bse = 0._wp
720     DO_2D_00_00
721
722        IF ( lconv(ji,jj) ) THEN
723           IF ( zdb_bl(ji,jj) > 0._wp ) THEN
724              zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
725              zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
726           ENDIF
727        ELSE
728           zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )
729           zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )
730        ENDIF
731     END_2D
732
733      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
734      !  Pycnocline gradients for scalars and velocity
735      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
736
737       DO_2D_00_00
738       !
739          IF ( lconv (ji,jj) ) THEN
740          ! Unstable conditions
741             IF( zdb_bl(ji,jj) > 0._wp ) THEN
742             ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero
743                ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) )
744                zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) )
745                zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) )
746                DO jk = 2 , ibld(ji,jj)
747                   znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
748                   zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
749                   zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
750                   zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
751                END DO
752             ENDIF
753          ELSE
754          ! stable conditions
755          ! if pycnocline profile only defined when depth steady of increasing.
756             IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady.
757                IF ( zdb_bl(ji,jj) > 0._wp ) THEN
758                  IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
759                      ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
760                      zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
761                      zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
762                      DO jk = 2, ibld(ji,jj)
763                         znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
764                         zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
765                         zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
766                         zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
767                      END DO
768                  ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
769                      ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
770                      zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
771                      zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
772                      DO jk = 2, ibld(ji,jj)
773                         znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
774                         zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
775                         zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
776                         zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
777                      END DO
778                   ENDIF ! IF (zhol >=0.5)
779                ENDIF    ! IF (zdb_bl> 0.)
780             ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero
781          ENDIF          ! IF (lconv)
782         !
783       END_2D
784!
785       DO_2D_00_00
786       !
787          IF ( lconv (ji,jj) ) THEN
788          ! Unstable conditions
789              zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / &
790            & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0)
791             zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / &
792           & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
793             DO jk = 2 , ibld(ji,jj)-1
794                znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
795                zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 )
796                zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
797             END DO
798          ELSE
799          ! stable conditions
800             zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
801             zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
802             DO jk = 2, ibld(ji,jj)
803                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
804                IF ( znd < 1.0 ) THEN
805                   zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
806                ELSE
807                   zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
808                ENDIF
809                zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
810             END DO
811          ENDIF
812         !
813       END_2D
814       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
815       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
816       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
817
818      ! WHERE ( lconv )
819      !     zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird
820      !     zvisml_sc = zdifml_sc
821      !     zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml )
822      !     zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml )
823      !     zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third
824      !     zbeta_v_sc = 1.0 -  2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml
825      !  ELSEWHERE
826      !     zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )
827      !     zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )
828      !  ENDWHERE
829       DO_2D_00_00
830          IF ( lconv(ji,jj) ) THEN
831            zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
832            zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
833            zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
834            zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
835            zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third
836            zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj)
837          ELSE
838            zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
839            zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
840         END IF
841       END_2D
842!
843       DO_2D_00_00
844          IF ( lconv(ji,jj) ) THEN
845             DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
846                 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
847                 !
848                 zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5
849                 !
850                 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) &
851                      &            *                                      ( 1.0 -               0.5 * zznd_ml**2 )
852             END DO
853             ! pycnocline - if present linear profile
854             IF ( zdh(ji,jj) > 0._wp ) THEN
855                DO jk = imld(ji,jj)+1 , ibld(ji,jj)
856                    zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
857                    !
858                    zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
859                    !
860                    zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
861                END DO
862             ENDIF
863             ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out
864             zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm)
865             ! could be taken out, take account of entrainment represents as a diffusivity
866             ! should remove w from here, represents entrainment
867          ELSE
868          ! stable conditions
869             DO jk = 2, ibld(ji,jj)
870                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
871                zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
872                zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
873             END DO
874          ENDIF   ! end if ( lconv )
875!
876       END_2D
877
878       !
879       ! calculate non-gradient components of the flux-gradient relationships
880       !
881! Stokes term in scalar flux, flux-gradient relationship
882       WHERE ( lconv )
883          zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln)
884          !
885          zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
886       ELSEWHERE
887          zsc_wth_1 = 2.0 * zwthav
888          !
889          zsc_ws_1 = 2.0 * zwsav
890       ENDWHERE
891
892
893       DO_2D_00_00
894         IF ( lconv(ji,jj) ) THEN
895           DO jk = 2, imld(ji,jj)
896              zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
897              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)
898              !
899              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)
900           END DO ! end jk loop
901         ELSE     ! else for if (lconv)
902 ! Stable conditions
903            DO jk = 2, ibld(ji,jj)
904               zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
905               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
906                    &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
907               !
908               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
909                    &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
910            END DO
911         ENDIF               ! endif for check on lconv
912
913       END_2D
914
915
916! 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)
917       WHERE ( lconv )
918          zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) )
919          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln )
920          zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )
921       ELSEWHERE
922          zsc_uw_1 = zustar**2
923          zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)
924       ENDWHERE
925
926       DO_2D_00_00
927          IF ( lconv(ji,jj) ) THEN
928             DO jk = 2, imld(ji,jj)
929                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
930                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   &
931                     &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) &
932                     &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) )
933!
934                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       &
935                     &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
936             END DO   ! end jk loop
937          ELSE
938! Stable conditions
939             DO jk = 2, ibld(ji,jj) ! corrected to ibld
940                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
941                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       &
942                     &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
943                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
944             END DO   ! end jk loop
945          ENDIF
946       END_2D
947
948! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
949
950       WHERE ( lconv )
951          zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
952          zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
953       ELSEWHERE
954          zsc_wth_1 = 0._wp
955          zsc_ws_1 = 0._wp
956       ENDWHERE
957
958       DO_2D_00_00
959          IF (lconv(ji,jj) ) THEN
960             DO jk = 2, imld(ji,jj)
961                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
962                ! calculate turbulent length scale
963                zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
964                     &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) )
965                zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
966                     &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
967                zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
968                ! non-gradient buoyancy terms
969                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 )
970                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 )
971             END DO
972          ELSE
973             DO jk = 2, ibld(ji,jj)
974                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
975                ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
976             END DO
977          ENDIF
978       END_2D
979
980
981       WHERE ( lconv )
982          zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
983          zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
984          zsc_vw_1 = 0._wp
985       ELSEWHERE
986         zsc_uw_1 = 0._wp
987         zsc_vw_1 = 0._wp
988       ENDWHERE
989
990       DO_2D_00_00
991          IF ( lconv(ji,jj) ) THEN
992             DO jk = 2 , imld(ji,jj)
993                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
994                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     &
995                     &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   &
996                     &                                          * zsc_uw_2(ji,jj)                                    )
997                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
998             END DO  ! jk loop
999          ELSE
1000          ! stable conditions
1001             DO jk = 2, ibld(ji,jj)
1002                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
1003                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
1004             END DO
1005          ENDIF
1006       END_2D
1007
1008! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
1009
1010       WHERE ( lconv )
1011          zsc_wth_1 = zwth0
1012          zsc_ws_1 = zws0
1013       ELSEWHERE
1014          zsc_wth_1 = 2.0 * zwthav
1015          zsc_ws_1 = zws0
1016       ENDWHERE
1017
1018       DO_2D_00_00
1019         IF ( lconv(ji,jj) ) THEN
1020            DO jk = 2, imld(ji,jj)
1021               zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
1022               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                &
1023                    &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
1024                    &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
1025                    &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) )
1026               !
1027               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  &
1028                    &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
1029                    &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
1030                    &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) )
1031            END DO
1032         ELSE
1033            DO jk = 2, ibld(ji,jj)
1034               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
1035               znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
1036               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
1037            &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
1038               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
1039            &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
1040            END DO
1041         ENDIF
1042       END_2D
1043
1044
1045       WHERE ( lconv )
1046          zsc_uw_1 = zustar**2
1047          zsc_vw_1 = ff_t * zustke * zhml
1048       ELSEWHERE
1049          zsc_uw_1 = zustar**2
1050          zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
1051          zsc_vw_1 = ff_t * zustke * zhbl
1052          zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1
1053       ENDWHERE
1054
1055       DO_2D_00_00
1056          IF ( lconv(ji,jj) ) THEN
1057            DO jk = 2, imld(ji,jj)
1058               zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
1059               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
1060               ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
1061                    & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)
1062               !
1063               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1064                    & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
1065            END DO
1066          ELSE
1067            DO jk = 2, ibld(ji,jj)
1068               znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
1069               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
1070               IF ( zznd_d <= 2.0 ) THEN
1071                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
1072                       &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
1073                  !
1074               ELSE
1075                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
1076                       & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
1077                  !
1078               ENDIF
1079
1080               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1081                    & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj)
1082               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1083                    & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
1084            END DO
1085          ENDIF
1086       END_2D
1087!
1088! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
1089
1090      DO_2D_00_00
1091         IF ( lconv(ji,jj) ) THEN
1092            DO jk = 2, ibld(ji,jj)
1093               znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
1094               IF ( znd >= 0.0 ) THEN
1095                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
1096                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
1097               ELSE
1098                  ghamu(ji,jj,jk) = 0._wp
1099                  ghamv(ji,jj,jk) = 0._wp
1100               ENDIF
1101            END DO
1102         ELSE
1103            DO jk = 2, ibld(ji,jj)
1104               znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
1105               IF ( znd >= 0.0 ) THEN
1106                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
1107                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
1108               ELSE
1109                  ghamu(ji,jj,jk) = 0._wp
1110                  ghamv(ji,jj,jk) = 0._wp
1111               ENDIF
1112            END DO
1113         ENDIF
1114      END_2D
1115
1116      ! pynocline contributions
1117       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small
1118       zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln )
1119       DO_2D_00_00
1120          DO jk= 2, ibld(ji,jj)
1121             znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
1122             ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
1123             ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
1124             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
1125             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)
1126             ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
1127          END DO
1128       END_2D
1129
1130! Entrainment contribution.
1131
1132       DO_2D_00_00
1133          IF ( lconv(ji,jj) ) THEN
1134            DO jk = 1, imld(ji,jj) - 1
1135               znd=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
1136               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
1137               ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
1138               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd
1139               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd
1140            END DO
1141            DO jk = imld(ji,jj), ibld(ji,jj)
1142               znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
1143               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
1144               ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
1145               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd )
1146               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd )
1147             END DO
1148          ENDIF
1149          ghamt(ji,jj,ibld(ji,jj)) = 0._wp
1150          ghams(ji,jj,ibld(ji,jj)) = 0._wp
1151          ghamu(ji,jj,ibld(ji,jj)) = 0._wp
1152          ghamv(ji,jj,ibld(ji,jj)) = 0._wp
1153       END_2D
1154
1155
1156       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1157       ! Need to put in code for contributions that are applied explicitly to
1158       ! the prognostic variables
1159       !  1. Entrainment flux
1160       !
1161       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1162
1163
1164
1165       ! rotate non-gradient velocity terms back to model reference frame
1166
1167       DO_2D_00_00
1168          DO jk = 2, ibld(ji,jj)
1169             ztemp = ghamu(ji,jj,jk)
1170             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
1171             ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
1172          END DO
1173       END_2D
1174
1175       IF(ln_dia_osm) THEN
1176          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1177       END IF
1178
1179! KPP-style Ri# mixing
1180       IF( ln_kpprimix) THEN
1181          DO_3D_10_10( 2, jpkm1 )
1182             z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   &
1183                  &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) &
1184                  &                 / (  e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) )
1185             z3dv(ji,jj,jk) = 0.5 * (  vv(ji,jj,jk-1,Kmm) -  vv(ji,jj  ,jk,Kmm) )   &
1186                  &                 * (  vv(ji,jj,jk-1,Kbb) -  vv(ji,jj  ,jk,Kbb) ) * wvmask(ji,jj,jk) &
1187                  &                 / (  e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) )
1188          END_3D
1189      !
1190         DO_3D_00_00( 2, jpkm1 )
1191            !                                          ! shear prod. at w-point weightened by mask
1192            zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
1193               &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
1194            !                                          ! local Richardson number
1195            zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
1196            zfri =  MIN( zri / rn_riinfty , 1.0_wp )
1197            zfri  = ( 1.0_wp - zfri * zfri )
1198            zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk)
1199         END_3D
1200
1201          DO_2D_00_00
1202             DO jk = ibld(ji,jj) + 1, jpkm1
1203                zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
1204                zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
1205             END DO
1206          END_2D
1207
1208       END IF ! ln_kpprimix = .true.
1209
1210! KPP-style set diffusivity large if unstable below BL
1211       IF( ln_convmix) THEN
1212          DO_2D_00_00
1213             DO jk = ibld(ji,jj) + 1, jpkm1
1214               IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
1215             END DO
1216          END_2D
1217       END IF ! ln_convmix = .true.
1218
1219       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
1220       CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1. )
1221
1222       ! GN 25/8: need to change tmask --> wmask
1223
1224     DO_3D_00_00( 2, jpkm1 )
1225          p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1226          p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
1227     END_3D
1228      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid  (sign unchanged), needed to caclulate gham[uv] on u and v grids
1229     CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1.,   &
1230      &                  ghamu, 'W', 1. , ghamv, 'W', 1. )
1231       DO_3D_00_00( 2, jpkm1 )
1232            ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
1233               &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1234
1235            ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
1236                &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1237
1238            ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk)
1239            ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk)
1240       END_3D
1241        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged)
1242        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged)
1243        CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1.,   &
1244         &                  ghamu, 'U', 1. , ghamv, 'V', 1. )
1245
1246       IF(ln_dia_osm) THEN
1247         SELECT CASE (nn_osm_wave)
1248         ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1).
1249         CASE(0:1)
1250            IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift
1251            IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift
1252            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke )
1253         ! Stokes drift read in from sbcwave  (=2).
1254         CASE(2)
1255            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift
1256            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift
1257            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* &
1258                 & SQRT(ut0sd**2 + vt0sd**2 ) )
1259         END SELECT
1260         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )            ! <Tw_NL>
1261         IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams )            ! <Sw_NL>
1262         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )            ! <uw_NL>
1263         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )            ! <vw_NL>
1264         IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 )            ! <Tw_0>
1265         IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0>
1266         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth
1267         IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli )               ! Initial boundary-layer depth
1268         IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth
1269         IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points
1270         IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc )         ! convective velocity scale
1271         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale
1272         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale
1273         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
1274         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke )
1275         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine
1276         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine
1277         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )               ! ML depth internal to zdf_osm routine
1278         IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine
1279         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )               ! ML depth internal to zdf_osm routine
1280         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )               ! ML depth internal to zdf_osm routine
1281         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )               ! average T in ML
1282      END IF
1283      ! Lateral boundary conditions on p_avt  (sign unchanged)
1284      CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1. )
1285      !
1286   END SUBROUTINE zdf_osm
1287
1288
1289   SUBROUTINE zdf_osm_init( Kmm ) 
1290     !!----------------------------------------------------------------------
1291     !!                  ***  ROUTINE zdf_osm_init  ***
1292     !!
1293     !! ** Purpose :   Initialization of the vertical eddy diffivity and
1294     !!      viscosity when using a osm turbulent closure scheme
1295     !!
1296     !! ** Method  :   Read the namosm namelist and check the parameters
1297     !!      called at the first timestep (nit000)
1298     !!
1299     !! ** input   :   Namlist namosm
1300     !!----------------------------------------------------------------------
1301     INTEGER, INTENT(in)    :: Kmm ! time level index (middle)
1302     !
1303     INTEGER  ::   ios            ! local integer
1304     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1305     !!
1306     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
1307          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
1308          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv
1309     !!----------------------------------------------------------------------
1310     !
1311     READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
1312901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
1313
1314     READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
1315902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
1316     IF(lwm) WRITE ( numond, namzdf_osm )
1317
1318     IF(lwp) THEN                    ! Control print
1319        WRITE(numout,*)
1320        WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
1321        WRITE(numout,*) '~~~~~~~~~~~~'
1322        WRITE(numout,*) '   Namelist namzdf_osm : set tke mixing parameters'
1323        WRITE(numout,*) '     Use namelist  rn_osm_la                     ln_use_osm_la = ', ln_use_osm_la
1324        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
1325        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
1326        WRITE(numout,*) '     Depth scale of Stokes drift                rn_osm_dstokes = ', rn_osm_dstokes
1327        WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
1328        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave
1329        SELECT CASE (nn_osm_wave)
1330        CASE(0)
1331           WRITE(numout,*) '     calculated assuming constant La#=0.3'
1332        CASE(1)
1333           WRITE(numout,*) '     calculated from Pierson Moskowitz wind-waves'
1334        CASE(2)
1335           WRITE(numout,*) '     calculated from ECMWF wave fields'
1336        END SELECT
1337        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
1338        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix
1339        WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
1340        WRITE(numout,*) '     maximum shear diffusivity at Rig = 0    (m2/s) rn_difri = ', rn_difri
1341        WRITE(numout,*) '     Use large mixing below BL when unstable       ln_convmix = ', ln_convmix
1342        WRITE(numout,*) '     diffusivity when unstable below BL     (m2/s) rn_difconv = ', rn_difconv
1343     ENDIF
1344
1345     !                              ! allocate zdfosm arrays
1346     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
1347
1348     call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl
1349
1350     IF( ln_zdfddm) THEN
1351        IF(lwp) THEN
1352           WRITE(numout,*)
1353           WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
1354           WRITE(numout,*) '    CAUTION : done in routine zdfosm, not in routine zdfddm '
1355        ENDIF
1356     ENDIF
1357
1358
1359     !set constants not in namelist
1360     !-----------------------------
1361
1362     IF(lwp) THEN
1363        WRITE(numout,*)
1364     ENDIF
1365
1366     IF (nn_osm_wave == 0) THEN
1367        dstokes(:,:) = rn_osm_dstokes
1368     END IF
1369
1370     ! Horizontal average : initialization of weighting arrays
1371     ! -------------------
1372
1373     SELECT CASE ( nn_ave )
1374
1375     CASE ( 0 )                ! no horizontal average
1376        IF(lwp) WRITE(numout,*) '          no horizontal average on avt'
1377        IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
1378        ! weighting mean arrays etmean
1379        !           ( 1  1 )
1380        ! avt = 1/4 ( 1  1 )
1381        !
1382        etmean(:,:,:) = 0.e0
1383
1384        DO_3D_00_00( 1, jpkm1 )
1385           etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
1386                &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
1387                &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
1388        END_3D
1389
1390     CASE ( 1 )                ! horizontal average
1391        IF(lwp) WRITE(numout,*) '          horizontal average on avt'
1392        ! weighting mean arrays etmean
1393        !           ( 1/2  1  1/2 )
1394        ! avt = 1/8 ( 1    2  1   )
1395        !           ( 1/2  1  1/2 )
1396        etmean(:,:,:) = 0.e0
1397
1398        DO_3D_00_00( 1, jpkm1 )
1399           etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
1400                & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
1401                &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
1402                &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
1403                &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
1404                &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
1405        END_3D
1406
1407     CASE DEFAULT
1408        WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
1409        CALL ctl_stop( ctmp1 )
1410
1411     END SELECT
1412
1413     ! Initialization of vertical eddy coef. to the background value
1414     ! -------------------------------------------------------------
1415     DO jk = 1, jpk
1416        avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1417     END DO
1418
1419     ! zero the surface flux for non local term and osm mixed layer depth
1420     ! ------------------------------------------------------------------
1421     ghamt(:,:,:) = 0.
1422     ghams(:,:,:) = 0.
1423     ghamu(:,:,:) = 0.
1424     ghamv(:,:,:) = 0.
1425     !
1426     IF( lwxios ) THEN
1427        CALL iom_set_rstw_var_active('wn')
1428        CALL iom_set_rstw_var_active('hbl')
1429        CALL iom_set_rstw_var_active('hbli')
1430     ENDIF
1431   END SUBROUTINE zdf_osm_init
1432
1433
1434   SUBROUTINE osm_rst( kt, Kmm, cdrw )
1435     !!---------------------------------------------------------------------
1436     !!                   ***  ROUTINE osm_rst  ***
1437     !!
1438     !! ** Purpose :   Read or write BL fields in restart file
1439     !!
1440     !! ** Method  :   use of IOM library. If the restart does not contain
1441     !!                required fields, they are recomputed from stratification
1442     !!----------------------------------------------------------------------
1443
1444     INTEGER         , INTENT(in) ::   kt     ! ocean time step index
1445     INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index (middle)
1446     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1447
1448     INTEGER ::   id1, id2   ! iom enquiry index
1449     INTEGER  ::   ji, jj, jk     ! dummy loop indices
1450     INTEGER  ::   iiki, ikt ! local integer
1451     REAL(wp) ::   zhbf           ! tempory scalars
1452     REAL(wp) ::   zN2_c           ! local scalar
1453     REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
1454     INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top)
1455     !!----------------------------------------------------------------------
1456     !
1457     !!-----------------------------------------------------------------------------
1458     ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
1459     !!-----------------------------------------------------------------------------
1460     IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
1461        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. )
1462        IF( id1 > 0 ) THEN                       ! 'wn' exists; read
1463           CALL iom_get( numror, jpdom_autoglo, 'wn', ww, ldxios = lrxios )
1464           WRITE(numout,*) ' ===>>>> :  ww read from restart file'
1465        ELSE
1466           ww(:,:,:) = 0._wp
1467           WRITE(numout,*) ' ===>>>> :  ww not in restart file, set to zero initially'
1468        END IF
1469        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
1470        id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. )
1471        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return
1472           CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
1473           CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli, ldxios = lrxios  )
1474           WRITE(numout,*) ' ===>>>> :  hbl & hbli read from restart file'
1475           RETURN
1476        ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate
1477           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
1478        END IF
1479     END IF
1480
1481     !!-----------------------------------------------------------------------------
1482     ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
1483     !!-----------------------------------------------------------------------------
1484     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
1485        IF(lwp) WRITE(numout,*) '---- osm-rst ----'
1486         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , ww  , ldxios = lwxios )
1487         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios )
1488         CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios )
1489        RETURN
1490     END IF
1491
1492     !!-----------------------------------------------------------------------------
1493     ! Getting hbl, no restart file with hbl, so calculate from surface stratification
1494     !!-----------------------------------------------------------------------------
1495     IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
1496     ALLOCATE( imld_rst(jpi,jpj) )
1497     ! w-level of the mixing and mixed layers
1498     CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )
1499     CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm)
1500     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
1501     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
1502     zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria
1503     !
1504     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
1505     DO_3D_11_11( 1, jpkm1 )
1506        ikt = mbkt(ji,jj)
1507        hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
1508        IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
1509     END_3D
1510     !
1511     DO_2D_11_11
1512        iiki = imld_rst(ji,jj)
1513        hbl (ji,jj) = gdepw(ji,jj,iiki  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth
1514     END_2D
1515     hbl = MAX(hbl,epsln)
1516     hbli(:,:) = hbl(:,:)
1517     DEALLOCATE( imld_rst )
1518     WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
1519   END SUBROUTINE osm_rst
1520
1521
1522   SUBROUTINE tra_osm( kt, Kmm, pts, Krhs )
1523      !!----------------------------------------------------------------------
1524      !!                  ***  ROUTINE tra_osm  ***
1525      !!
1526      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
1527      !!
1528      !! ** Method  :   ???
1529      !!----------------------------------------------------------------------
1530      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
1531      !!----------------------------------------------------------------------
1532      INTEGER                                  , INTENT(in)    :: kt        ! time step index
1533      INTEGER                                  , INTENT(in)    :: Kmm, Krhs ! time level indices
1534      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation
1535      !
1536      INTEGER :: ji, jj, jk
1537      !
1538      IF( kt == nit000 ) THEN
1539         IF(lwp) WRITE(numout,*)
1540         IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
1541         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1542      ENDIF
1543
1544      IF( l_trdtra )   THEN                    !* Save ta and sa trends
1545         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
1546         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
1547      ENDIF
1548
1549      ! add non-local temperature and salinity flux
1550      DO_3D_00_00( 1, jpkm1 )
1551         pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      &
1552            &                 - (  ghamt(ji,jj,jk  )  &
1553            &                    - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm)
1554         pts(ji,jj,jk,jp_sal,Krhs) =  pts(ji,jj,jk,jp_sal,Krhs)                      &
1555            &                 - (  ghams(ji,jj,jk  )  &
1556            &                    - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
1557      END_3D
1558
1559
1560      ! save the non-local tracer flux trends for diagnostic
1561      IF( l_trdtra )   THEN
1562         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
1563         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
1564!!bug gm jpttdzdf ==> jpttosm
1565         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt )
1566         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds )
1567         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
1568      ENDIF
1569
1570      IF(sn_cfctl%l_prtctl) THEN
1571         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm  - Ta: ', mask1=tmask,   &
1572         &             tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
1573      ENDIF
1574      !
1575   END SUBROUTINE tra_osm
1576
1577
1578   SUBROUTINE trc_osm( kt )          ! Dummy routine
1579      !!----------------------------------------------------------------------
1580      !!                  ***  ROUTINE trc_osm  ***
1581      !!
1582      !! ** Purpose :   compute and add to the passive tracer trend the non-local
1583      !!                 passive tracer flux
1584      !!
1585      !!
1586      !! ** Method  :   ???
1587      !!----------------------------------------------------------------------
1588      !
1589      !!----------------------------------------------------------------------
1590      INTEGER, INTENT(in) :: kt
1591      WRITE(*,*) 'trc_osm: Not written yet', kt
1592   END SUBROUTINE trc_osm
1593
1594
1595   SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs )
1596      !!----------------------------------------------------------------------
1597      !!                  ***  ROUTINE dyn_osm  ***
1598      !!
1599      !! ** Purpose :   compute and add to the velocity trend the non-local flux
1600      !! copied/modified from tra_osm
1601      !!
1602      !! ** Method  :   ???
1603      !!----------------------------------------------------------------------
1604      INTEGER                             , INTENT( in )  ::  kt          ! ocean time step index
1605      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
1606      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
1607      !
1608      INTEGER :: ji, jj, jk   ! dummy loop indices
1609      !!----------------------------------------------------------------------
1610      !
1611      IF( kt == nit000 ) THEN
1612         IF(lwp) WRITE(numout,*)
1613         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
1614         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1615      ENDIF
1616      !code saving tracer trends removed, replace with trdmxl_oce
1617
1618      DO_3D_00_00( 1, jpkm1 )
1619         puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs)                      &
1620            &                 - (  ghamu(ji,jj,jk  )  &
1621            &                    - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm)
1622         pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs)                      &
1623            &                 - (  ghamv(ji,jj,jk  )  &
1624            &                    - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm)
1625      END_3D
1626      !
1627      ! code for saving tracer trends removed
1628      !
1629   END SUBROUTINE dyn_osm
1630
1631   !!======================================================================
1632END MODULE zdfosm
Note: See TracBrowser for help on using the repository browser.