source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfosm.F90 @ 10030

Last change on this file since 10030 was 10030, checked in by gm, 3 years ago

#1911 (ENHANCE-04): RK3 branch - step II.3 remove e3uw_$ e3vw_$, except e3.w_0 and use only after e3 in dyn/trazdf

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