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

Last change on this file since 13237 was 13237, checked in by smasson, 3 months ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

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