source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfosm.F90 @ 12616

Last change on this file since 12616 was 12616, checked in by techene, 11 months ago

all: add e3 substitute (sometimes it requires to add ze3t/u/v/w) and limit precompiled files lines to about 130 character, OCE/ASM/asminc.F90, OCE/DOM/domzgr_substitute.h90, OCE/ISF/isfcpl.F90, OCE/SBC/sbcice_cice, OCE/CRS/crsini.F90 : add key_LF

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