source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90 @ 9106

Last change on this file since 9106 was 9106, checked in by cetlod, 3 years ago

Bugfix on ZDF

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