source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfosm.F90 @ 13517

Last change on this file since 13517 was 13517, checked in by hadcv, 5 months ago

Tiling for modules after tra_ldf

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