source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/ZDF/zdfosm.F90 @ 10297

Last change on this file since 10297 was 10297, checked in by smasson, 23 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of mppmin/max/sum, see #2133

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