New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfkpp.F90 in branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 10 years ago

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169

  • Property svn:keywords set to Id
File size: 79.5 KB
Line 
1MODULE zdfkpp
2   !!======================================================================
3   !!                       ***  MODULE  zdfkpp  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the KPP
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !! History :  OPA  ! 2000-03 (W.G. Large, J. Chanut) Original code
8   !!            8.1  ! 2002-06 (J.M. Molines) for real case CLIPPER 
9   !!            8.2  ! 2003-10 (Chanut J.) re-writting
10   !!   NEMO     1.0  ! 2005-01 (C. Ethe, G. Madec) Free form, F90 + creation of tra_kpp routine
11   !!            3.3  ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
12   !!----------------------------------------------------------------------
13#if defined key_zdfkpp   ||   defined key_esopa
14   !!----------------------------------------------------------------------
15   !!   'key_zdfkpp'                                             KPP scheme
16   !!----------------------------------------------------------------------
17   !!   zdf_kpp       : update momentum and tracer Kz from a kpp scheme
18   !!   zdf_kpp_init  : initialization, namelist read, and parameters control
19   !!   tra_kpp       : compute and add to the T & S trend the non-local flux
20   !!   trc_kpp       : compute and add to the passive tracer trend the non-local flux (lk_top=T)
21   !!----------------------------------------------------------------------
22   USE oce            ! ocean dynamics and active tracers
23   USE dom_oce        ! ocean space and time domain
24   USE zdf_oce        ! ocean vertical physics
25   USE sbc_oce        ! surface boundary condition: ocean
26   USE phycst         ! physical constants
27   USE eosbn2         ! equation of state
28   USE zdfddm         ! double diffusion mixing
29   USE in_out_manager ! I/O manager
30   USE lib_mpp        ! MPP library
31   USE wrk_nemo       ! work arrays
32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
33   USE prtctl         ! Print control
34   USE trdmod_oce     ! ocean trends definition
35   USE trdtra         ! tracers trends
36   USE timing         ! Timing
37   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   zdf_kpp       ! routine called by step.F90
43   PUBLIC   zdf_kpp_init  ! routine called by opa.F90
44   PUBLIC   tra_kpp       ! routine called by step.F90
45   PUBLIC   trc_kpp       ! routine called by trcstp.F90
46
47   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfkpp = .TRUE.    !: KPP vertical mixing flag
48
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghats    !: non-local scalar mixing term (gamma/<ws>o)
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wt0      !: surface temperature flux for non local flux
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ws0      !: surface salinity flux for non local flux
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hkpp     !: boundary layer depth
53
54   !                          !!* Namelist namzdf_kpp *
55   REAL(wp) ::   rn_difmiw     ! constant internal wave viscosity (m2/s)
56   REAL(wp) ::   rn_difsiw     ! constant internal wave diffusivity (m2/s)
57   REAL(wp) ::   rn_riinfty    ! local Richardson Number limit for shear instability
58   REAL(wp) ::   rn_difri      ! maximum shear mixing at Rig = 0    (m2/s)
59   REAL(wp) ::   rn_bvsqcon    ! Brunt-Vaisala squared (1/s**2) for maximum convection
60   REAL(wp) ::   rn_difcon     ! maximum mixing in interior convection (m2/s)
61   INTEGER  ::   nn_ave        ! = 0/1 flag for horizontal average on avt, avmu, avmv
62
63#if defined key_zdfddm
64   !                                        !!! ** Double diffusion Mixing
65   REAL(wp) ::   difssf  = 1.e-03_wp         ! maximum salt fingering mixing
66   REAL(wp) ::   Rrho0   = 1.9_wp            ! limit for salt  fingering mixing
67   REAL(wp) ::   difsdc  = 1.5e-06_wp        ! maximum diffusive convection mixing
68#endif
69   LOGICAL  ::   ln_kpprimix       ! Shear instability mixing
70
71   !                                        !!! ** General constants  **
72   REAL(wp) ::   epsln   = 1.0e-20_wp        ! a small positive number   
73   REAL(wp) ::   pthird  = 1._wp/3._wp       ! 1/3
74   REAL(wp) ::   pfourth = 1._wp/4._wp       ! 1/4
75
76   !                                        !!! ** Boundary Layer Turbulence Parameters  **
77   REAL(wp) ::   vonk     = 0.4_wp           ! von Karman's constant
78   REAL(wp) ::   epsilon  = 0.1_wp           ! nondimensional extent of the surface layer
79   REAL(wp) ::   rconc1   = 5.0_wp           ! standard flux profile function parmaeters
80   REAL(wp) ::   rconc2   = 16.0_wp          !         "        "
81   REAL(wp) ::   rconcm   = 8.38_wp          ! momentum flux profile fit
82   REAL(wp) ::   rconam   = 1.26_wp          !         "       "
83   REAL(wp) ::   rzetam   = -.20_wp          !         "       "       
84   REAL(wp) ::   rconcs   = 98.96_wp         !  scalar  flux profile fit
85   REAL(wp) ::   rconas   = -28.86_wp        !         "       "
86   REAL(wp) ::   rzetas   = -1.0_wp          !         "       " 
87   
88   !                                        !!! ** Boundary Layer Depth Diagnostic  **
89   REAL(wp) ::   Ricr     = 0.3_wp           ! critical bulk Richardson Number
90   REAL(wp) ::   rcekman  = 0.7_wp           ! coefficient for ekman depth 
91   REAL(wp) ::   rcmonob  = 1.0_wp           ! coefficient for Monin-Obukhov depth
92   REAL(wp) ::   rconcv   = 1.7_wp           ! ratio of interior buoyancy frequency to its value at entrainment depth
93   REAL(wp) ::   hbf      = 1.0_wp           ! fraction of bound. layer depth to which absorbed solar
94      !                                      ! rad. and contributes to surf. buo. forcing
95   REAL(wp) ::   Vtc                         ! function of rconcv,rconcs,epsilon,vonk,Ricr
96   
97   !                                        !!! ** Nonlocal Boundary Layer Mixing **
98   REAL(wp) ::   rcstar   = 5.0_wp           ! coefficient for convective nonlocal transport
99   REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s   
100   REAL(wp) ::   rcg                         ! non-dimensional coefficient for nonlocal transport
101
102#if ! defined key_kppcustom
103   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   del     ! array for reference mean values of vertical integration
104#endif
105
106#if defined key_kpplktb
107   !                                         !!! ** Parameters for lookup table for turbulent velocity scales **
108   INTEGER, PARAMETER ::   nilktb   = 892     ! number of values for zehat in KPP lookup table
109   INTEGER, PARAMETER ::   njlktb   = 482     ! number of values for ustar in KPP lookup table
110   INTEGER, PARAMETER ::   nilktbm1 = nilktb-1   !
111   INTEGER, PARAMETER ::   njlktbm1 = njlktb-1   !
112
113   REAL(wp), DIMENSION(nilktb,njlktb) ::   wmlktb   ! lookup table for the turbulent vertical velocity scale (momentum)
114   REAL(wp), DIMENSION(nilktb,njlktb) ::   wslktb   ! lookup table for the turbulent vertical velocity scale (tracers)
115
116   REAL(wp) ::   dehatmin = -4.e-7_wp    ! minimum limit for zhat in lookup table (m3/s3)
117   REAL(wp) ::   dehatmax = 0._wp        ! maximum limit for zhat in lookup table (m3/s3)
118   REAL(wp) ::   ustmin   = 0._wp        ! minimum limit for ustar in lookup table (m/s)
119   REAL(wp) ::   ustmax   = 0.04_wp      ! maximum limit for ustar in lookup table (m/s)   
120   REAL(wp) ::   dezehat                 ! delta zhat in lookup table
121   REAL(wp) ::   deustar                 ! delta ustar in lookup table
122#endif
123   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   ratt   ! attenuation coef  (already defines in module traqsr,
124   !                                    ! but only if the solar radiation penetration is considered)
125   
126   !                                    !!! * penetrative solar radiation coefficient *
127   REAL(wp) ::   rabs = 0.58_wp          ! fraction associated with xsi1
128   REAL(wp) ::   xsi1 = 0.35_wp          ! first depth of extinction
129   REAL(wp) ::   xsi2 = 23.0_wp          ! second depth of extinction
130      !                           ! (default values: water type Ib)
131
132   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean, eumean, evmean   ! coeff. used for hor. smoothing at t-, u- & v-points
133       
134#if defined key_c1d
135   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rig    !: gradient Richardson number
136   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rib    !: bulk Richardson number
137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   buof   !: buoyancy forcing
138   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mols   !: moning-Obukhov length scale
139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ekdp   !: Ekman depth
140#endif
141
142   INTEGER  ::   jip = 62 , jjp = 111
143
144   !! * Substitutions
145#  include "domzgr_substitute.h90"
146#  include "vectopt_loop_substitute.h90"
147#  include  "zdfddm_substitute.h90"
148   !!----------------------------------------------------------------------
149   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
150   !! $Id$
151   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
152   !!----------------------------------------------------------------------
153CONTAINS
154
155   INTEGER FUNCTION zdf_kpp_alloc()
156      !!----------------------------------------------------------------------
157      !!                 ***  FUNCTION zdf_kpp_alloc  ***
158      !!----------------------------------------------------------------------
159      ALLOCATE( ghats(jpi,jpj,jpk), wt0(jpi,jpj), ws0(jpi,jpj), hkpp(jpi,jpj), &
160#if ! defined key_kpplktb
161         &      del(jpk,jpk),                                                  &
162#endif
163         &      ratt(jpk),                                                     &
164         &      etmean(jpi,jpj,jpk), eumean(jpi,jpj,jpk), evmean(jpi,jpj,jpk), &
165#if defined key_c1d
166         &      rig (jpi,jpj,jpk), rib(jpi,jpj,jpk), buof(jpi,jpj,jpk),        &
167         &      mols(jpi,jpj,jpk), ekdp(jpi,jpj),                              &
168#endif
169         &      STAT= zdf_kpp_alloc )
170         !
171      IF( lk_mpp             )   CALL mpp_sum ( zdf_kpp_alloc )
172      IF( zdf_kpp_alloc /= 0 )   CALL ctl_warn('zdf_kpp_alloc: failed to allocate arrays')
173   END FUNCTION zdf_kpp_alloc
174
175
176   SUBROUTINE zdf_kpp( kt )
177      !!----------------------------------------------------------------------
178      !!                   ***  ROUTINE zdf_kpp  ***
179      !!
180      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
181      !!      coefficients and non local mixing using K-profile parameterization
182      !!
183      !! ** Method :   The boundary layer depth hkpp is diagnosed at tracer points
184      !!      from profiles of buoyancy, and shear, and the surface forcing.
185      !!      Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
186      !!
187      !!                      Kx =  hkpp  Wx(sigma) G(sigma) 
188      !!
189      !!             and the non local term ghat = Cs / Ws(sigma) / hkpp
190      !!      Below hkpp  the coefficients are the sum of mixing due to internal waves
191      !!      shear instability and double diffusion.
192      !!
193      !!      -1- Compute the now interior vertical mixing coefficients at all depths.
194      !!      -2- Diagnose the boundary layer depth.
195      !!      -3- Compute the now boundary layer vertical mixing coefficients.
196      !!      -4- Compute the now vertical eddy vicosity and diffusivity.
197      !!      -5- Smoothing
198      !!
199      !!        N.B. The computation is done from jk=2 to jpkm1
200      !!             Surface value of avt avmu avmv are set once a time to zero
201      !!             in routine zdf_kpp_init.
202      !!
203      !! ** Action  :   update the non-local terms ghats
204      !!                update avt, avmu, avmv (before vertical eddy coef.)
205      !!
206      !! References : Large W.G., Mc Williams J.C. and Doney S.C.             
207      !!         Reviews of Geophysics, 32, 4, November 1994
208      !!         Comments in the code refer to this paper, particularly
209      !!         the equation number. (LMD94, here after)
210      !!----------------------------------------------------------------------
211      USE oce     , zviscos => ua   ! temp. array for viscosities use ua as workspace
212      USE oce     , zdiffut => va   ! temp. array for diffusivities use sa as workspace
213      !!
214      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
215      !!
216      INTEGER ::   ji, jj, jk              ! dummy loop indices
217      INTEGER ::   ikbot, jkmax, jkm1, jkp2   !
218
219      REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube     !
220      REAL(wp) ::   zrhos, zalbet, zbeta, zthermal, zhalin, zatt1   !
221      REAL(wp) ::   zref, zt, zs, zh, zu, zv, zrh                   ! Bulk richardson number
222      REAL(wp) ::   zrib, zrinum, zdVsq, zVtsq                      !
223      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm   ! Velocity scales
224#if defined key_kpplktb
225      INTEGER ::    il, jl                                          ! Lookup table or Analytical functions
226      REAL(wp) ::   ud, zfrac, ufrac, zwam, zwbm, zwas, zwbs        !
227#else
228      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      !
229#endif
230      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density
231#if ! defined key_kppcustom     
232      INTEGER  ::   jm                          ! dummy loop indices
233      REAL(wp) ::   zr1, zr2, zr3, zr4, zrhop   ! Compression terms
234#endif
235      REAL(wp) ::   zflag, ztemp, zrn2, zdep21, zdep32, zdep43
236      REAL(wp) ::   zdku2, zdkv2, ze3sqr, zsh2, zri, zfri          ! Interior richardson mixing
237      REAL(wp), POINTER, DIMENSION(:,:) ::   zmoek                 ! Moning-Obukov limitation
238      REAL(wp), POINTER, DIMENSION(:)   ::   zmoa, zekman               
239      REAL(wp)                          ::   zmob, zek
240      REAL(wp), POINTER, DIMENSION(:,:) ::   zdepw, zdift, zvisc   ! The pipe
241      REAL(wp), POINTER, DIMENSION(:,:) ::   zdept
242      REAL(wp), POINTER, DIMENSION(:,:) ::   zriblk
243      REAL(wp), POINTER, DIMENSION(:)   ::   zhmax, zria, zhbl 
244      REAL(wp) ::   zflagri, zflagek, zflagmo, zflagh, zflagkb   !
245      REAL(wp), POINTER, DIMENSION(:)   ::   za2m, za3m, zkmpm, za2t, za3t, zkmpt   ! Shape function (G)
246      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
247#if defined key_zdfddm
248      REAL(wp) ::   zrrau, zds, zavdds, zavddt,zinr   ! double diffusion mixing
249      REAL(wp), POINTER, DIMENSION(:,:)   ::     zdifs
250      REAL(wp), POINTER, DIMENSION(:)     ::   za2s, za3s, zkmps
251      REAL(wp) ::                            zkm1s
252      REAL(wp), POINTER, DIMENSION(:,:)   ::   zblcs
253      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdiffus
254#endif
255      REAL(wp), POINTER, DIMENSION(:,:) ::   zBo, zBosol, zustar         ! Surface buoyancy forcing, friction velocity
256      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask, zblcm, zblct
257      !!--------------------------------------------------------------------
258      !
259      IF( nn_timing == 1 )  CALL timing_start('zdf_kpp')
260      !
261      CALL wrk_alloc( jpi, zmoa, zekman, zhmax, zria, zhbl )
262      CALL wrk_alloc( jpi, za2m, za3m, zkmpm, za2t, za3t, zkmpt )
263      CALL wrk_alloc( jpi,2, zriblk )
264      CALL wrk_alloc( jpi,3, zmoek, kjstart = 0 )
265      CALL wrk_alloc( jpi,3, zdept )
266      CALL wrk_alloc( jpi,4, zdepw, zdift, zvisc )
267      CALL wrk_alloc( jpi,jpj, zBo, zBosol, zustar )
268      CALL wrk_alloc( jpi,jpk, zmask, zblcm, zblct )
269#if defined key_zdfddm
270      CALL wrk_alloc( jpi,4, zdifs )
271      CALL wrk_alloc( jpi, zmoa, za2s, za3s, zkmps )
272      CALL wrk_alloc( jpi,jpk, zblcs )
273      CALL wrk_alloc( jpi,jpi,jpk, zdiffus )
274#endif
275
276      zviscos(:,:,:) = 0.
277      zblcm  (:,:  ) = 0. 
278      zdiffut(:,:,:) = 0.
279      zblct  (:,:  ) = 0. 
280#if defined key_zdfddm
281      zdiffus(:,:,:) = 0.
282      zblcs  (:,:  ) = 0. 
283#endif
284      ghats(:,:,:) = 0.
285     
286      zBo   (:,:) = 0.
287      zBosol(:,:) = 0.
288      zustar(:,:) = 0.
289
290
291      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
292      ! I. Interior diffusivity and viscosity at w points ( T interfaces)
293      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
294      DO jk = 2, jpkm1
295         DO jj = 2, jpjm1
296            DO ji = fs_2, fs_jpim1 
297               ! Mixing due to internal waves breaking
298               ! -------------------------------------
299               avmu(ji,jj,jk)  = rn_difmiw 
300               avt (ji,jj,jk)  = rn_difsiw             
301               ! Mixing due to vertical shear instability
302               ! -------------------------------------               
303               IF( ln_kpprimix ) THEN         
304                  ! Compute the gradient Richardson  number at interfaces (zri):
305                  ! LMD94, eq. 27 (is vertical smoothing needed : Rig=N^2 / (dz(u))^2 + (dz(v))^2
306                  zdku2 =   ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
307                     &    * ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
308                     &    + ( un(ji    ,jj,jk - 1) - un(ji    ,jj,jk) ) &
309                     &    * ( un(ji    ,jj,jk - 1) - un(ji    ,jj,jk) ) 
310                 
311                  zdkv2 =   ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
312                     &    * ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
313                     &    + ( vn(ji,    jj,jk - 1) - vn(ji,    jj,jk) ) &
314                     &    * ( vn(ji,    jj,jk - 1) - vn(ji,    jj,jk) ) 
315
316                  ze3sqr = 1. / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
317                  ! Square of vertical shear  at interfaces
318                  zsh2   = 0.5 * ( zdku2 + zdkv2 ) * ze3sqr
319                  zri    = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + epsln ) 
320#if defined key_c1d
321                  ! save the gradient richardson number
322                  rig(ji,jj,jk) = zri * tmask(ji,jj,jk)
323#endif                 
324                  ! Evaluate f of Ri (zri) for shear instability store in zfri
325                  ! LMD94, eq. 28a,b,c, figure 3 ; Rem: p1 is 3, hard coded
326                  zfri  = MAX( zri , 0. )
327                  zfri  = MIN( zfri / rn_riinfty , 1.0 )
328                  zfri  = ( 1.0 - zfri * zfri )
329                  zfri  = zfri * zfri  * zfri
330                  ! add shear contribution to mixing coef.
331                  avmu(ji,jj,jk) =  avmu(ji,jj,jk) + rn_difri * zfri   
332                  avt (ji,jj,jk) =  avt (ji,jj,jk) + rn_difri * zfri   
333               ENDIF
334#if defined key_zdfddm 
335               avs (ji,jj,jk) =  avt (ji,jj,jk)             
336               !  Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90
337               ! ------------------------------------------------------------------
338               ! only retains positive value of rrau
339               zrrau = MAX( rrau(ji,jj,jk), epsln )
340               zds   = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal)
341               IF( zrrau > 1. .AND. zds > 0.) THEN
342                  !
343                  ! Salt fingering case.
344                  !---------------------
345                  ! Compute interior diffusivity for double diffusive mixing of
346                  ! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
347                  ! After that set interior diffusivity for double diffusive mixing
348                  ! of temperature
349                  zavdds = MIN( zrrau, Rrho0 )
350                  zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 )
351                  zavdds = 1.0 - zavdds * zavdds 
352                  zavdds = zavdds * zavdds * zavdds 
353                  zavdds = difssf * zavdds 
354                  zavddt = 0.7 * zavdds
355               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN
356                  !
357                  ! Diffusive convection case.
358                  !---------------------------
359                  ! Compute interior diffusivity for double diffusive mixing of
360                  ! temperature (Marmorino and Caldwell, 1976);
361                  ! Compute interior diffusivity for double diffusive mixing of salinity
362                  zinr   = 1. / zrrau
363                  zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) ) 
364                  zavddt = difsdc * zavddt
365                  IF( zrrau < 0.5) THEN
366                     zavdds = zavddt * 0.15 * zrrau
367                  ELSE
368                     zavdds = zavddt * (1.85 * zrrau - 0.85 ) 
369                  ENDIF
370               ELSE
371                  zavddt = 0.
372                  zavdds = 0.
373               ENDIF 
374               ! Add double diffusion contribution to temperature and salinity  mixing coefficients.
375               avt (ji,jj,jk) =  avt (ji,jj,jk) +  zavddt 
376               avs (ji,jj,jk) =  avs (ji,jj,jk) +  zavdds         
377#endif                     
378            END DO
379         END DO
380      END DO
381
382
383      ! Radiative (zBosol) and non radiative (zBo) surface buoyancy
384      !JMM at the time zdfkpp is called, q still holds the sum q + qsr
385      !---------------------------------------------------------------------
386      DO jj = 2, jpjm1
387         DO ji = fs_2, fs_jpim1     
388            IF( nn_eos < 1) THEN   
389               zt     = tsn(ji,jj,1,jp_tem)
390               zs     = tsn(ji,jj,1,jp_sal) - 35.0
391               zh     = fsdept(ji,jj,1)
392               !  potential volumic mass
393               zrhos  = rhop(ji,jj,1)
394               zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta
395                  &                               - 0.203814e-03 ) * zt   &
396                  &                               + 0.170907e-01 ) * zt   &
397                  &   + 0.665157e-01                                      &
398                  &   +     ( - 0.678662e-05 * zs                         &
399                  &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
400                  &   +   ( ( - 0.302285e-13 * zh                         &
401                  &           - 0.251520e-11 * zs                         &
402                  &           + 0.512857e-12 * zt * zt           ) * zh   &
403                  &           - 0.164759e-06 * zs                         &
404                  &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
405                  &                               + 0.380374e-04 ) * zh
406
407               zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta
408                  &                            - 0.301985e-05 ) * zt      &
409                  &   + 0.785567e-03                                      &
410                  &   + (     0.515032e-08 * zs                           &
411                  &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
412                  &   +(  (   0.121551e-17 * zh                           &
413                  &         - 0.602281e-15 * zs                           &
414                  &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
415                  &                             + 0.408195e-10   * zs     &
416                  &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
417                  &                             - 0.121555e-07 ) * zh
418
419               zthermal = zbeta * zalbet / ( rcp * zrhos + epsln )
420               zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs
421            ELSE
422               zrhos    = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) )
423               zthermal = rn_alpha / ( rcp * zrhos + epsln )
424               zhalin   = rn_beta * tsn(ji,jj,1,jp_sal) * rcs
425               zbeta    = rn_beta
426            ENDIF
427            ! Radiative surface buoyancy force
428            zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
429            ! Non radiative surface buoyancy force
430            zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * ( emp(ji,jj)-rnf(ji,jj) )  &
431               &                                         -  grav * zbeta * rcs * sfx(ji,jj)
432            ! Surface Temperature flux for non-local term
433            wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* r1_rau0_rcp * tmask(ji,jj,1)
434            ! Surface salinity flux for non-local term
435            ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)                          &
436               &             + sfx(ji,jj)                                     ) * rcs * tmask(ji,jj,1) 
437         ENDDO
438      ENDDO
439
440      zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. ) 
441      !  Compute surface buoyancy forcing, Monin Obukhov and Ekman depths 
442      !------------------------------------------------------------------   
443      DO jj = 2, jpjm1
444         DO ji = fs_2, fs_jpim1
445            !  Reference surface density = density at first T point level   
446            zrhos         = rhop(ji,jj,1) + zflageos * rau0 * ( 1. - tmask(ji,jj,1) ) 
447            ! Friction velocity (zustar), at T-point : LMD94 eq. 2
448            zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos +  epsln ) )
449         ENDDO
450      ENDDO
451
452!CDIR NOVERRCHK 
453      !                                               ! ===============
454      DO jj = 2, jpjm1                                 !  Vertical slab
455         !                                             ! ===============
456         
457         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
458         ! II Compute Boundary layer mixing coef. and diagnose the new boundary layer depth
459         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
460         
461         ! Initialization
462         jkmax       = 0
463         zdept (:,:) = 0.
464         zdepw (:,:) = 0.
465         zriblk(:,:) = 0.
466         zmoek (:,:) = 0.
467         zvisc (:,:) = 0.
468         zdift (:,:) = 0.
469#if defined key_zdfddm
470         zdifs (:,:) = 0.
471#endif
472         zmask (:,:) = 0.
473         DO ji = fs_2, fs_jpim1
474            zria(ji ) = 0.
475            ! Maximum boundary layer depth
476            ikbot     = mbkt(ji,jj)    ! ikbot is the last T point in the water
477            zhmax(ji) = fsdept(ji,jj,ikbot) - 0.001     
478            ! Compute Monin obukhov length scale at the surface and Ekman depth:
479            zbuofdep   = zBo(ji,jj) + zBosol(ji,jj) * ratt(1)
480            zekman(ji) = rcekman * zustar(ji,jj) / ( ABS( ff(ji,jj) ) + epsln )
481            zucube     = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) 
482            zmoa(ji)   = zucube / ( vonk * ( zbuofdep + epsln ) )   
483#if defined key_c1d
484            ! store the surface buoyancy forcing
485            zstabl        = 0.5 + SIGN( 0.5, zbuofdep )
486            buof(ji,jj,1) = zbuofdep * tmask(ji,jj,1)
487            ! store the moning-oboukov length scale at surface
488            zmob          = zstabl * zmoa(ji) + ( 1.0 - zstabl ) * fsdept(ji,jj,1)
489            mols(ji,jj,1) = MIN( zmob , zhmax(ji) ) * tmask(ji,jj,1)
490            ! store Ekman depth
491            zek           = zstabl * zekman(ji) + ( 1.0 - zstabl ) * fsdept(ji,jj,1) 
492            ekdp(ji,jj )  = MIN( zek , zhmax(ji) ) * tmask(ji,jj,1) 
493#endif
494         END DO     
495         ! Compute the pipe
496         ! ---------------------
497         DO jk = 2, jpkm1
498            DO ji = fs_2, fs_jpim1
499               ! Compute bfsfc = Bo + radiative contribution down to hbf*depht
500               zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * ratt(jk)
501               ! Flag (zstabl  = 1) if positive forcing
502               zstabl   =  0.5 + SIGN(  0.5, zbuofdep)
503
504               !   Compute bulk richardson number zrib at depht
505               !-------------------------------------------------------
506               !                           [Br - B(d)] * d         zrinum
507               !             Rib(z) = ----------------------- = -------------
508               !                       |Vr - V(d)|^2 + Vt(d)^2   zdVsq + zVtsq
509               !
510               ! First compute zt,zs,zu,zv = means in the surface layer < epsilon*depht 
511               ! Else surface values are taken at the first T level.
512               ! For stability, resolved vertical shear is computed with "before velocities".
513               zref = epsilon * fsdept(ji,jj,jk)
514#if defined key_kppcustom
515               ! zref = gdept(1)
516               zref = fsdept(ji,jj,1)
517               zt   = tsn(ji,jj,1,jp_tem)
518               zs   = tsn(ji,jj,1,jp_sal)
519               zrh  = rhop(ji,jj,1)
520               zu   = ( ub(ji,jj,1) + ub(ji - 1,jj    ,1) ) / MAX( 1. , umask(ji,jj,1) + umask(ji - 1,jj   ,1) )
521               zv   = ( vb(ji,jj,1) + vb(ji    ,jj - 1,1) ) / MAX( 1. , vmask(ji,jj,1) + vmask(ji   ,jj - 1,1) )
522#else
523               zt   = 0.
524               zs   = 0.
525               zu   = 0.
526               zv   = 0.
527               zrh  = 0.
528               ! vertically integration over the upper epsilon*gdept(jk) ; del () array is computed once in zdf_kpp_init
529               DO jm = 1, jpkm1
530                  zt   = zt  + del(jk,jm) * tsn(ji,jj,jm,jp_tem)
531                  zs   = zs  + del(jk,jm) * tsn(ji,jj,jm,jp_sal)
532                  zu   = zu  + 0.5 * del(jk,jm) &
533                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
534                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
535                  zv   = zv  + 0.5 * del(jk,jm) &
536                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
537                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
538                  zrh  = zrh + del(jk,jm) * rhop(ji,jj,jm)
539               END DO
540#endif
541               zsr = SQRT( ABS( tsn(ji,jj,jk,jp_sal) ) )
542               ! depth
543               zh = fsdept(ji,jj,jk)
544               ! compute compression terms on density
545               ze  = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6
546               zbw = (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4
547               zb  = zbw + ze * zs
548               
549               zd  = -2.042967e-2
550               zc  =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896
551               zaw = ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788
552               za  = ( zd*zsr + zc ) *zs + zaw
553               
554               zb1 =   (-0.1909078*zt+7.390729 ) *zt-55.87545
555               za1 = ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077
556               zkw = ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6
557               zk0 = ( zb1*zsr + za1 )*zs + zkw
558               zcomp =   1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )
559               
560#if defined key_kppcustom
561               ! potential density of water(zrh = zt,zs at level jk):
562               zrhdr = zrh / zcomp
563#else
564               ! potential density of water(ztref,zsref at level jk):
565               ! compute volumic mass pure water at atm pressure
566               IF ( nn_eos < 1 ) THEN
567                  zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt   &
568                     &       -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594
569                  ! seawater volumic mass atm pressure
570                  zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt   &
571                     &   -4.0899e-3 ) *zt+0.824493
572                  zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3
573                  zr4= 4.8314e-4             
574                  ! potential volumic mass (reference to the surface)
575                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1                 
576                  zrhdr = zrhop / zcomp
577               ELSE
578                  zrhdr = zrh / zcomp
579               ENDIF
580#endif
581               
582               ! potential density of ambiant water at level jk :
583               zrhd   = ( rhd(ji,jj,jk) * rau0 + rau0 ) 
584               
585               ! And now the Rib number numerator .
586               zrinum = grav * ( zrhd - zrhdr ) / rau0
587               zrinum = zrinum * ( fsdept(ji,jj,jk) - zref ) * tmask(ji,jj,jk)
588           
589               ! Resolved shear contribution to Rib at depth T-point (zdVsq)
590               ztx    =   ( ub( ji , jj ,jk)   +  ub(ji - 1, jj ,jk) ) &
591                  &     / MAX( 1. , umask( ji , jj ,jk) + umask(ji - 1, jj ,jk) )   
592               zty    =   ( vb( ji , jj ,jk)   +  vb(ji  ,jj - 1,jk) ) &
593                  &     / MAX( 1., vmask( ji , jj ,jk) + vmask(ji  ,jj - 1,jk) ) 
594               
595               zdVsq  = ( zu - ztx ) * ( zu - ztx ) + ( zv - zty ) * ( zv - zty )
596               
597               ! Scalar turbulent velocity scale zws for hbl=gdept
598               zscale = zstabl + ( 1.0 - zstabl ) * epsilon
599               zehat  = vonk * zscale * fsdept(ji,jj,jk) * zbuofdep
600               zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)             
601               zeta   = zehat / ( zucube + epsln )
602               
603               IF( zehat > 0. ) THEN
604                  ! Stable case
605                  zws  = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
606               ELSE
607                  ! Unstable case
608#if defined key_kpplktb
609                  ! use lookup table
610                  zd     = zehat - dehatmin
611                  il     = INT( zd / dezehat )
612                  il     = MIN( il, nilktbm1 )
613                  il     = MAX( il, 1 )
614                 
615                  ud     = zustar(ji,jj) - ustmin
616                  jl     = INT( ud / deustar )
617                  jl     = MIN( jl, njlktbm1 )
618                  jl     = MAX( jl, 1 )
619                 
620                  zfrac  = zd / dezehat - FLOAT( il ) 
621                  ufrac  = ud / deustar - FLOAT( jl )
622                  zwas   = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
623                  zwbs   = ( 1. - zfrac ) * wslktb(il,jl  ) + zfrac * wslktb(il+1,jl  )
624                  !
625                  zws    = ( 1. - ufrac ) * zwbs + ufrac * zwas
626#else
627                  ! use analytical functions:
628                  zcons  = 0.5 + SIGN( 0.5 , ( rzetas - zeta ) )
629                  zwcons = vonk * zustar(ji,jj) * ( ( ABS( rconas - rconcs * zeta ) )**pthird ) 
630                  zwsun  = vonk * zustar(ji,jj) * SQRT( ABS ( 1.0 - rconc2 * zeta ) )
631                  !
632                  zws    = zcons * zwcons +  ( 1.0 - zcons) * zwsun
633#endif
634               ENDIF
635               
636               ! Turbulent shear contribution to Rib (zVtsq) bv frequency at levels  ( ie T-point jk)
637               zrn2   = 0.5 * ( rn2(ji,jj,jk) + rn2(ji,jj,jk+1) )   
638               zbvzed = SQRT( ABS( zrn2 ) ) 
639               zVtsq  = fsdept(ji,jj,jk) * zws * zbvzed  * Vtc
640               
641               ! Finally, the bulk Richardson number at depth fsdept(i,j,k)
642               zrib  = zrinum   / ( zdVsq + zVtsq + epsln )
643 
644               ! Find subscripts around the boundary layer depth, build the pipe
645               ! ----------------------------------------------------------------
646
647               ! Flag (zflagri = 1) if zrib < Ricr 
648               zflagri = 0.5 + SIGN( 0.5, ( Ricr - zrib ) )
649               !  Flag (zflagh  = 1) if still within overall boundary layer
650               zflagh  = 0.5 + SIGN( 0.5, ( fsdept(ji,jj,1) - zdept(ji,2) ) )
651               
652               ! Ekman layer depth
653               zek     = zstabl * zekman(ji) + ( 1.0 - zstabl ) * zhmax(ji)
654               zflag   = 0.5 + SIGN( 0.5, ( zek - fsdept(ji,jj,jk-1) ) )
655               zek     = zflag * zek + ( 1.0 - zflag ) * zhmax(ji)
656               zflagek = 0.5 + SIGN( 0.5, ( zek - fsdept(ji,jj,jk) ) )
657               ! Flag (zflagmo = 1) if still within stable Monin-Obukhov and in water
658               zmob    = zucube / ( vonk * ( zbuofdep + epsln ) ) 
659               ztemp   = zstabl * zmob + ( 1.0 - zstabl) * zhmax(ji)
660               ztemp   = MIN( ztemp , zhmax(ji) ) 
661               zflagmo = 0.5 + SIGN( 0.5, ( ztemp - fsdept(ji,jj,jk) ) )             
662
663               ! No limitation by Monin Obukhov or Ekman depths:
664!               zflagek = 1.0
665!               zflagmo = 0.5 + SIGN( 0.5, ( zhmax(ji) - fsdept(ji,jj,jk) ) )
666
667               ! Load  pipe via zflagkb  for later calculations
668               ! Flag (zflagkb = 1) if zflagh = 1 and (zflagri = 0 or zflagek = 0 or zflagmo = 0)
669               zflagkb = zflagh * ( 1.0 - ( zflagri * zflagek * zflagmo ) )
670
671               zmask(ji,jk) = zflagh
672               jkp2         = MIN( jk+2 , ikbot )
673               jkm1         = MAX( jk-1 , 2 )
674               jkmax        = MAX( jkmax, jk * INT( REAL( zflagh+epsln ) ) ) 
675
676               zdept(ji,1)  = zdept(ji,1) + zflagkb * fsdept(ji,jj,jk-1) 
677               zdept(ji,2)  = zdept(ji,2) + zflagkb * fsdept(ji,jj,jk  ) 
678               zdept(ji,3)  = zdept(ji,3) + zflagkb * fsdept(ji,jj,jk+1) 
679
680               zdepw(ji,1)  = zdepw(ji,1) + zflagkb * fsdepw(ji,jj,jk-1) 
681               zdepw(ji,2)  = zdepw(ji,2) + zflagkb * fsdepw(ji,jj,jk  ) 
682               zdepw(ji,3)  = zdepw(ji,3) + zflagkb * fsdepw(ji,jj,jk+1)
683               zdepw(ji,4)  = zdepw(ji,4) + zflagkb * fsdepw(ji,jj,jkp2) 
684
685               zriblk(ji,1) = zriblk(ji,1) + zflagkb * zria(ji)
686               zriblk(ji,2) = zriblk(ji,2) + zflagkb * zrib
687
688               zmoek (ji,0) = zmoek (ji,0) + zflagkb * zek
689               zmoek (ji,1) = zmoek (ji,1) + zflagkb * zmoa(ji)
690               zmoek (ji,2) = zmoek (ji,2) + zflagkb * ztemp 
691               ! Save Monin Obukhov depth
692               zmoa  (ji)   = zmob
693           
694               zvisc(ji,1) = zvisc(ji,1) + zflagkb * avmu(ji,jj,jkm1)
695               zvisc(ji,2) = zvisc(ji,2) + zflagkb * avmu(ji,jj,jk  )
696               zvisc(ji,3) = zvisc(ji,3) + zflagkb * avmu(ji,jj,jk+1)
697               zvisc(ji,4) = zvisc(ji,4) + zflagkb * avmu(ji,jj,jkp2)
698               
699               zdift(ji,1) = zdift(ji,1) + zflagkb * avt (ji,jj,jkm1)
700               zdift(ji,2) = zdift(ji,2) + zflagkb * avt (ji,jj,jk  )
701               zdift(ji,3) = zdift(ji,3) + zflagkb * avt (ji,jj,jk+1)
702               zdift(ji,4) = zdift(ji,4) + zflagkb * avt (ji,jj,jkp2)
703
704#if defined key_zdfddm
705               zdifs(ji,1) = zdifs(ji,1) + zflagkb * avs (ji,jj,jkm1)
706               zdifs(ji,2) = zdifs(ji,2) + zflagkb * avs (ji,jj,jk  )
707               zdifs(ji,3) = zdifs(ji,3) + zflagkb * avs (ji,jj,jk+1)
708               zdifs(ji,4) = zdifs(ji,4) + zflagkb * avs (ji,jj,jkp2)
709#endif               
710               ! Save the Richardson number
711               zria  (ji)   = zrib 
712#if defined key_c1d
713               ! store buoyancy length scale
714               buof(ji,jj,jk) = zbuofdep * tmask(ji,jj,jk) 
715               ! store Monin Obukhov
716               zmob           = zstabl * zmob + ( 1.0 - zstabl) * fsdept(ji,jj,1)
717               mols(ji,jj,jk) = MIN( zmob , zhmax(ji) ) * tmask(ji,jj,jk) 
718               ! Bulk Richardson number
719               rib(ji,jj,jk)  = zrib * tmask(ji,jj,jk)             
720#endif               
721            END DO
722         END DO
723         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
724         ! III PROCESS THE PIPE
725         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
726         
727         DO ji = fs_2, fs_jpim1 
728           
729            ! Find the boundary layer depth zhbl
730            ! ----------------------------------------
731           
732            ! Interpolate monin Obukhov and critical Ri mumber depths   
733            ztemp = zdept(ji,2) - zdept(ji,1)
734            zflag = ( Ricr - zriblk(ji,1) ) / ( zriblk(ji,2) - zriblk(ji,1)  + epsln )
735            zhrib = zdept(ji,1) + zflag * ztemp     
736
737            IF( zriblk(ji,2) < Ricr ) zhrib = zhmax(ji) 
738         
739            IF( zmoek(ji,2) < zdept(ji,2) ) THEN
740               IF ( zmoek(ji,1) < 0. ) THEN
741                  zmob = zdept(ji,2) - epsln
742               ELSE
743                  zmob = ztemp + zmoek(ji,1) - zmoek(ji,2)
744                  zmob = ( zmoek(ji,1) * zdept(ji,2) - zmoek(ji,2) * zdept(ji,1) ) / zmob
745                  zmob = MAX( zmob , zdept(ji,1) + epsln )               
746               ENDIF
747            ELSE           
748               zmob = zhmax(ji) 
749            ENDIF
750            ztemp   = MIN( zmob , zmoek(ji,0) )
751                         
752            ! Finally, the boundary layer depth, zhbl
753            zhbl(ji) = MAX( fsdept(ji,jj,1) + epsln, MIN( zhrib , ztemp ) )
754           
755            ! Save hkpp for further diagnostics (optional)
756            hkpp(ji,jj) = zhbl(ji) * tmask(ji,jj,1) 
757         
758            ! Correct mask if zhbl < fsdepw(ji,jj,2) for no viscosity/diffusivity enhancement at fsdepw(ji,jj,2)
759            !     zflag = 1 if zhbl(ji) > fsdepw(ji,jj,2)
760            IF( zhbl(ji) < fsdepw(ji,jj,2) ) zmask(ji,2) = 0.
761         
762           
763            !  Velocity scales at depth zhbl
764            ! -----------------------------------
765           
766            !  Compute bouyancy forcing down to zhbl
767            ztemp    = -hbf * zhbl(ji)
768            zatt1    = 1.0 - ( rabs * EXP( ztemp / xsi1 ) + ( 1.0 - rabs ) * EXP( ztemp / xsi2 ) )
769            zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * zatt1
770            zstabl   = 0.5 + SIGN( 0.5 , zbuofdep ) 
771
772            zbuofdep = zbuofdep + zstabl * epsln
773
774            zscale = zstabl + ( 1.0 - zstabl ) * epsilon         
775            zehat  = vonk * zscale * zhbl(ji) * zbuofdep
776            zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)             
777            zeta   = zehat / ( zucube + epsln )
778           
779            IF( zehat > 0. ) THEN
780               ! Stable case
781               zws  = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
782               zwm  = zws
783            ELSE
784               ! Unstable case
785#if defined key_kpplktb
786               ! use lookup table
787               zd     = zehat - dehatmin
788               il     = INT( zd / dezehat )
789               il     = MIN( il, nilktbm1 )
790               il     = MAX( il, 1 )
791               
792               ud     = zustar(ji,jj) - ustmin
793               jl     = INT( ud / deustar )
794               jl     = MIN( jl, njlktbm1 )
795               jl     = MAX( jl, 1 )
796               
797               zfrac  = zd / dezehat - FLOAT( il ) 
798               ufrac  = ud / deustar - FLOAT( jl )
799               zwas   = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
800               zwbs   = ( 1. - zfrac ) * wslktb(il,jl  ) + zfrac * wslktb(il+1,jl  )
801               zwam   = ( 1. - zfrac ) * wmlktb(il,jl+1) + zfrac * wmlktb(il+1,jl+1)
802               zwbm   = ( 1. - zfrac ) * wmlktb(il,jl  ) + zfrac * wmlktb(il+1,jl  )
803               !
804               zws    = ( 1. - ufrac ) * zwbs + ufrac * zwas
805               zwm    = ( 1. - ufrac ) * zwbm + ufrac * zwam
806#else
807               ! use analytical functions
808               zconm  = 0.5 + SIGN( 0.5, ( rzetam - zeta) )
809               zcons  = 0.5 + SIGN( 0.5, ( rzetas - zeta) )
810               
811               ! Momentum : zeta < rzetam (zconm = 1)
812               ! Scalars  : zeta < rzetas (zcons = 1)
813               zwconm = zustar(ji,jj) * vonk * ( ( ABS( rconam - rconcm * zeta) )**pthird )
814               zwcons = zustar(ji,jj) * vonk * ( ( ABS( rconas - rconcs * zeta) )**pthird )
815               
816               ! Momentum : rzetam <= zeta < 0 (zconm = 0)
817               ! Scalars  : rzetas <= zeta < 0 (zcons = 0) 
818               zwmun  = SQRT( ABS( 1.0 - rconc2 * zeta ) )
819               zwsun  = vonk * zustar(ji,jj) * zwmun
820               zwmun  = vonk * zustar(ji,jj) * SQRT(zwmun)
821               !
822               zwm    = zconm * zwconm + ( 1.0 - zconm ) * zwmun
823               zws    = zcons * zwcons + ( 1.0 - zcons ) * zwsun
824               
825#endif
826            ENDIF
827           
828           
829            ! Viscosity, diffusivity values and derivatives at h
830            ! --------------------------------------------------------
831           
832            ! check between at which interfaces is located zhbl(ji)
833            ! ztemp = 1, zdepw(ji,2) < zhbl <  zdepw(ji,3)
834            ! ztemp = 0, zdepw(ji,1) < zhbl <  zdepw(ji,2)
835            ztemp  =  0.5 + SIGN( 0.5, ( zhbl(ji) - zdepw(ji,2) ) ) 
836            zdep21 =   zdepw(ji,2) - zdepw(ji,1) + epsln
837            zdep32 =   zdepw(ji,3) - zdepw(ji,2) + epsln
838            zdep43 =   zdepw(ji,4) - zdepw(ji,3) + epsln 
839           
840            ! Compute R as in LMD94, eq D5b
841            zdelta =  ( zhbl(ji) - zdepw(ji,2) ) *         ztemp   / zdep32   &
842               &    + ( zhbl(ji) - zdepw(ji,1) ) * ( 1.0 - ztemp ) / zdep21 
843           
844            ! Compute the vertical derivative of viscosities (zdzh) at z=zhbl(ji)
845            zdzup  =  ( zvisc(ji,2) - zvisc(ji,3) ) *         ztemp   / zdep32 &
846               &    + ( zvisc(ji,1) - zvisc(ji,2) ) * ( 1.0 - ztemp ) / zdep21
847           
848            zdzdn  =  ( zvisc(ji,3) - zvisc(ji,4) ) *         ztemp   / zdep43 &
849               &    + ( zvisc(ji,2) - zvisc(ji,3) ) * ( 1.0 - ztemp ) / zdep32
850           
851            ! LMD94, eq D5b :         
852            zdzh   = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
853            zdzh   = MAX( zdzh , 0. )
854           
855            ! Compute viscosities (zvath) at z=zhbl(ji), LMD94 eq D5a
856            zvath  =          ztemp   * ( zvisc(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
857               &    + ( 1.0 - ztemp ) * ( zvisc(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
858           
859            ! Compute G (zgat1) and its derivative (zdat1) at z=hbl(ji), LMD94 eq 18
860           
861            ! Vertical derivative of velocity scale divided by velocity scale squared at z=hbl(ji)
862            ! (non zero only in stable conditions)
863            zflag  =  -zstabl * rconc1 * zbuofdep / ( zucube * zustar(ji,jj) + epsln )
864           
865            ! G at its derivative at z=hbl:
866            zgat1  = zvath  / ( zhbl(ji) * ( zwm + epsln )  )
867            zdat1  = -zdzh  / ( zwm + epsln ) -  zflag * zvath / zhbl(ji)
868           
869            ! G coefficients, LMD94 eq 17
870            za2m(ji) = -2.0 + 3.0 * zgat1 - zdat1
871            za3m(ji) =  1.0 - 2.0 * zgat1 + zdat1
872
873           
874            ! Compute the vertical derivative of temperature diffusivities (zdzh) at z=zhbl(ji)
875            zdzup  =  ( zdift(ji,2) - zdift(ji,3) ) *         ztemp   / zdep32 &
876               &    + ( zdift(ji,1) - zdift(ji,2) ) * ( 1.0 - ztemp ) / zdep21
877           
878            zdzdn  =  ( zdift(ji,3) - zdift(ji,4) ) *         ztemp   / zdep43 &
879               &    + ( zdift(ji,2) - zdift(ji,3) ) * ( 1.0 - ztemp ) / zdep32
880           
881            ! LMD94, eq D5b :         
882            zdzh   = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
883            zdzh   = MAX( zdzh , 0. )
884           
885           
886            ! Compute diffusivities (zvath) at z=zhbl(ji), LMD94 eq D5a
887            zvath  =          ztemp   * ( zdift(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
888               &    + ( 1.0 - ztemp ) * ( zdift(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
889                       
890            ! G at its derivative at z=hbl:
891            zgat1  = zvath  / ( zhbl(ji) * ( zws + epsln )  )
892            zdat1  = -zdzh  / ( zws + epsln ) -  zflag * zvath / zhbl(ji)
893           
894            ! G coefficients, LMD94 eq 17
895            za2t(ji) = -2.0 + 3.0 * zgat1 - zdat1
896            za3t(ji) =  1.0 - 2.0 * zgat1 + zdat1
897
898#if defined key_zdfddm
899            ! Compute the vertical derivative of salinities diffusivities (zdzh) at z=zhbl(ji)
900            zdzup  =  ( zdifs(ji,2) - zdifs(ji,3) ) *         ztemp   / zdep32 &
901               &    + ( zdifs(ji,1) - zdifs(ji,2) ) * ( 1.0 - ztemp ) / zdep21
902           
903            zdzdn  =  ( zdifs(ji,3) - zdifs(ji,4) ) *         ztemp   / zdep43 &
904               &    + ( zdifs(ji,2) - zdifs(ji,3) ) * ( 1.0 - ztemp ) / zdep32
905           
906            ! LMD94, eq D5b :         
907            zdzh   = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
908            zdzh   = MAX( zdzh , 0. )           
909           
910            ! Compute diffusivities (zvath) at z=zhbl(ji), LMD94 eq D5a
911            zvath  =          ztemp   * ( zdifs(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
912               &    + ( 1.0 - ztemp ) * ( zdifs(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
913                       
914            ! G at its derivative at z=hbl:
915            zgat1  = zvath  / ( zhbl(ji) * ( zws + epsln )  )
916            zdat1  = -zdzh  / ( zws + epsln ) -  zflag * zvath / zhbl(ji)
917           
918            ! G coefficients, LMD94 eq 17
919            za2s(ji) = -2.0 + 3.0 * zgat1 - zdat1
920            za3s(ji) =  1.0 - 2.0 * zgat1 + zdat1
921#endif
922
923            !-------------------turn off interior matching here------
924            !          za2(ji,1) = -2.0
925            !          za3(ji,1) =  1.0
926            !          za2(ji,2) = -2.0
927            !          za3(ji,2) =  1.0
928            !--------------------------------------------------------
929           
930            !  Compute Enhanced Mixing Coefficients (LMD94,eq D6)
931            ! ---------------------------------------------------------------
932           
933            ! Delta
934            zdelta  = ( zhbl(ji)  - zdept(ji,1) ) / ( zdept(ji,2) - zdept(ji,1) + epsln )
935            zdelta2 = zdelta * zdelta
936           
937            !  Mixing coefficients at first level above h (zdept(ji,1))
938            ! and at first interface in the pipe (zdepw(ji,2))
939           
940            ! At first T level above h (zdept(ji,1)) (always in the boundary layer)
941            zsig    = zdept(ji,1) / zhbl(ji)
942            ztemp   = zstabl * zsig  + ( 1.0 - zstabl ) * MIN( zsig , epsilon )
943            zehat   = vonk * ztemp * zhbl(ji) * zbuofdep
944            zeta    = zehat / ( zucube + epsln)
945            zwst    = vonk * zustar(ji,jj) / ( ABS( 1.0 + rconc1 * zeta ) + epsln)
946            zwm     = zstabl * zwst + ( 1.0 - zstabl ) * zwm
947            zws     = zstabl * zwst + ( 1.0 - zstabl ) * zws
948
949            zkm1m  = zhbl(ji) * zwm * zsig * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
950            zkm1t  = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
951#if defined key_zdfddm
952            zkm1s  = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
953#endif                       
954            ! At first W level in the pipe (zdepw(ji,2)) (not always in the boundary layer ):
955            zsig    = MIN( zdepw(ji,2) / zhbl(ji) , 1.0 )
956            ztemp   = zstabl * zsig + ( 1.0 - zstabl ) * MIN( zsig , epsilon )
957            zehat   = vonk * ztemp * zhbl(ji) * zbuofdep
958            zeta    = zehat / ( zucube + epsln )
959            zwst    = vonk * zustar(ji,jj) / ( ABS( 1.0 + rconc1 * zeta ) + epsln)
960            zws     = zstabl * zws + ( 1.0 - zstabl ) * zws
961            zwm     = zstabl * zws + ( 1.0 - zstabl ) * zwm
962
963            zkmpm(ji) = zhbl(ji) * zwm * zsig * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
964            zkmpt(ji) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
965#if defined key_zdfddm
966            zkmps(ji) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
967#endif 
968     
969            ! check if this point is in the boundary layer,else take interior viscosity/diffusivity:
970            zflag       = 0.5 + SIGN( 0.5, ( zhbl(ji) - zdepw(ji,2) ) )
971            zkmpm(ji) = zkmpm(ji) * zflag + ( 1.0 - zflag ) * zvisc(ji,2)
972            zkmpt(ji) = zkmpt(ji) * zflag + ( 1.0 - zflag ) * zdift(ji,2)
973#if defined key_zdfddm
974            zkmps(ji) = zkmps(ji) * zflag + ( 1.0 - zflag ) * zdifs(ji,2)
975#endif
976
977            ! Enhanced viscosity/diffusivity at zdepw(ji,2)
978            ztemp     = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1m + zdelta2 * zkmpm(ji)
979            zkmpm(ji) = ( 1.0 - zdelta ) * zvisc(ji,2) + zdelta * ztemp
980            ztemp     = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1t + zdelta2 * zkmpt(ji)
981            zkmpt(ji) = ( 1.0 - zdelta ) * zdift(ji,2) + zdelta * ztemp
982#if defined key_zdfddm
983            ztemp     = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1s + zdelta2 * zkmps(ji)
984            zkmps(ji) = ( 1.0 - zdelta ) * zdifs(ji,2) + zdelta * ztemp
985#endif           
986
987         END DO
988         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
989         ! IV. Compute vertical eddy viscosity and diffusivity coefficients
990         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
991         
992         DO jk  = 2, jkmax
993           
994            ! Compute turbulent velocity scales on the interfaces
995            ! --------------------------------------------------------
996            DO  ji = fs_2, fs_jpim1
997               zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * zatt1
998               zstabl   = 0.5 + SIGN( 0.5 , zbuofdep ) 
999               zbuofdep = zbuofdep + zstabl * epsln         
1000               zsig    = fsdepw(ji,jj,jk) / zhbl(ji)
1001               ztemp   = zstabl * zsig + ( 1. - zstabl ) * MIN( zsig , epsilon )
1002               zehat   = vonk * ztemp * zhbl(ji) * zbuofdep
1003               zucube  = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)
1004               zeta    = zehat / ( zucube + epsln )
1005
1006               IF( zehat > 0. ) THEN
1007                  ! Stable case
1008                  zws  = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
1009                  zwm  = zws
1010               ELSE
1011                  ! Unstable case
1012#if defined key_kpplktb
1013                  ! use lookup table
1014                  zd     = zehat - dehatmin
1015                  il     = INT( zd / dezehat )
1016                  il     = MIN( il, nilktbm1 )
1017                  il     = MAX( il, 1 )
1018                 
1019                  ud     = zustar(ji,jj) - ustmin
1020                  jl     = INT( ud / deustar )
1021                  jl     = MIN( jl, njlktbm1 )
1022                  jl     = MAX( jl, 1 )
1023                 
1024                  zfrac  = zd / dezehat - FLOAT( il ) 
1025                  ufrac  = ud / deustar - FLOAT( jl )
1026                  zwas   = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
1027                  zwbs   = ( 1. - zfrac ) * wslktb(il,jl  ) + zfrac * wslktb(il+1,jl  )
1028                  zwam   = ( 1. - zfrac ) * wmlktb(il,jl+1) + zfrac * wmlktb(il+1,jl+1)
1029                  zwbm   = ( 1. - zfrac ) * wmlktb(il,jl  ) + zfrac * wmlktb(il+1,jl  )
1030                  !
1031                  zws    = ( 1. - ufrac ) * zwbs + ufrac * zwas
1032                  zwm    = ( 1. - ufrac ) * zwbm + ufrac * zwam
1033#else
1034                  ! use analytical functions
1035                  zconm  = 0.5 + SIGN( 0.5, ( rzetam - zeta) )
1036                  zcons  = 0.5 + SIGN( 0.5, ( rzetas - zeta) )
1037                 
1038                  ! Momentum : zeta < rzetam (zconm = 1)
1039                  ! Scalars  : zeta < rzetas (zcons = 1)
1040                  zwconm = zustar(ji,jj) * vonk * ( ( ABS( rconam - rconcm * zeta) )**pthird )
1041                  zwcons = zustar(ji,jj) * vonk * ( ( ABS( rconas - rconcs * zeta) )**pthird )
1042                 
1043                  ! Momentum : rzetam <= zeta < 0 (zconm = 0)
1044                  ! Scalars  : rzetas <= zeta < 0 (zcons = 0) 
1045                  zwmun  = SQRT( ABS( 1.0 - rconc2 * zeta ) )
1046                  zwsun  = vonk * zustar(ji,jj) * zwmun
1047                  zwmun  = vonk * zustar(ji,jj) * SQRT(zwmun)
1048                  !
1049                  zwm    = zconm * zwconm + ( 1.0 - zconm ) * zwmun
1050                  zws    = zcons * zwcons + ( 1.0 - zcons ) * zwsun
1051                 
1052#endif
1053               ENDIF
1054               
1055               zblcm(ji,jk) = zhbl(ji) * zwm * zsig  * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
1056               zblct(ji,jk) = zhbl(ji) * zws * zsig  * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
1057#if defined key_zdfddm
1058               zblcs(ji,jk) = zhbl(ji) * zws * zsig  * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
1059#endif             
1060               !  Compute Nonlocal transport term = ghats * <ws>o
1061               ! ----------------------------------------------------
1062               ghats(ji,jj,jk-1) = ( 1. - zstabl ) * rcg / ( zws * zhbl(ji) + epsln ) * tmask(ji,jj,jk)
1063
1064            END DO
1065         END DO     
1066         !  Combine interior and boundary layer coefficients and nonlocal term
1067         ! -----------------------------------------------------------------------
1068         DO jk = 2, jpkm1   
1069            DO ji = fs_2, fs_jpim1
1070               zflag = zmask(ji,jk) * zmask(ji,jk+1)
1071               zviscos(ji,jj,jk) = ( 1.0 - zmask(ji,jk) )         * avmu (ji,jj,jk) & ! interior viscosities
1072                  &              +                        zflag   * zblcm(ji,jk    ) & ! boundary layer viscosities
1073                  &              + zmask(ji,jk) * ( 1.0 - zflag ) * zkmpm(ji       )   ! viscosity enhancement at W_level near zhbl
1074               
1075               zviscos(ji,jj,jk) = zviscos(ji,jj,jk) * tmask(ji,jj,jk)   
1076
1077           
1078               zdiffut(ji,jj,jk) = ( 1.0 - zmask(ji,jk) )          * avt (ji,jj,jk) & ! interior diffusivities
1079                  &              +                        zflag   * zblct(ji,jk   ) & ! boundary layer diffusivities
1080                  &              + zmask(ji,jk) * ( 1.0 - zflag ) * zkmpt(ji      )   ! diffusivity enhancement at W_level near zhbl
1081                       
1082               zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) * tmask(ji,jj,jk) 
1083#if defined key_zdfddm
1084               zdiffus(ji,jj,jk) = ( 1.0 - zmask(ji,jk) )         * avs (ji,jj,jk) & ! interior diffusivities
1085                  &              +                        zflag   * zblcs(ji,jk   ) & ! boundary layer diffusivities
1086                  &              + zmask(ji,jk) * ( 1.0 - zflag ) * zkmps(ji      )   ! diffusivity enhancement at W_level near zhbl
1087               zdiffus(ji,jj,jk) = zdiffus(ji,jj,jk) * tmask(ji,jj,jk) 
1088#endif               
1089               ! Non local flux in the boundary layer only
1090               ghats(ji,jj,jk-1) = zmask(ji,jk) * ghats(ji,jj,jk-1)
1091
1092            ENDDO
1093         END DO
1094         !                                             ! ===============
1095      END DO                                           !   End of slab
1096      !                                                ! ===============
1097
1098      ! Lateral boundary conditions on zvicos and zdiffus  (sign unchanged)
1099      CALL lbc_lnk( zviscos(:,:,:), 'U', 1. )  ; CALL lbc_lnk( zdiffut(:,:,:), 'W', 1. ) 
1100#if defined key_zdfddm 
1101      CALL lbc_lnk( zdiffus(:,:,:), 'W', 1. ) 
1102#endif
1103
1104      SELECT CASE ( nn_ave )
1105         !
1106         CASE ( 0 )             ! no viscosity and diffusivity smoothing
1107
1108            DO jk = 2, jpkm1
1109               DO jj = 2, jpjm1
1110                  DO ji = fs_2, fs_jpim1
1111                     avmu(ji,jj,jk) = ( zviscos(ji,jj,jk) + zviscos(ji+1,jj,jk) ) &
1112                        &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1113                     
1114                     avmv(ji,jj,jk) = ( zviscos(ji,jj,jk) + zviscos(ji,jj+1,jk) ) &
1115                        &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1116                     
1117                     avt (ji,jj,jk) =  zdiffut(ji,jj,jk) * tmask(ji,jj,jk) 
1118#if defined key_zdfddm     
1119                     avs (ji,jj,jk) =  zdiffus(ji,jj,jk) * tmask(ji,jj,jk) 
1120#endif
1121                  END DO
1122               END DO
1123            END DO
1124           
1125         CASE ( 1 )                ! viscosity and diffusivity smoothing
1126            !                     
1127            !           ( 1/2  1  1/2 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
1128            ! avt = 1/8 ( 1    2  1   )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
1129            !           ( 1/2  1  1/2 )              ( 1/2  1/2 )
1130 
1131            DO jk = 2, jpkm1
1132               DO jj = 2, jpjm1
1133                  DO ji = fs_2, fs_jpim1
1134
1135                     avmu(ji,jj,jk) = (      zviscos(ji  ,jj  ,jk) + zviscos(ji+1,jj  ,jk)   &
1136                        &              +.5*( zviscos(ji  ,jj-1,jk) + zviscos(ji+1,jj-1,jk)   &
1137                        &                   +zviscos(ji  ,jj+1,jk) + zviscos(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
1138                     
1139                     avmv(ji,jj,jk) = (      zviscos(ji  ,jj  ,jk) + zviscos(ji  ,jj+1,jk)   &
1140                        &              +.5*( zviscos(ji-1,jj  ,jk) + zviscos(ji-1,jj+1,jk)   &
1141                        &                   +zviscos(ji+1,jj  ,jk) + zviscos(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
1142 
1143                     avt (ji,jj,jk) = ( .5*( zdiffut(ji-1,jj+1,jk) + zdiffut(ji-1,jj-1,jk)    &
1144                        &                   +zdiffut(ji+1,jj+1,jk) + zdiffut(ji+1,jj-1,jk) )  &
1145                        &              +1.*( zdiffut(ji-1,jj  ,jk) + zdiffut(ji  ,jj+1,jk)    &
1146                        &                   +zdiffut(ji  ,jj-1,jk) + zdiffut(ji+1,jj  ,jk) )  &
1147                        &              +2.*  zdiffut(ji  ,jj  ,jk)                          ) * etmean(ji,jj,jk)
1148#if defined key_zdfddm   
1149                     avs (ji,jj,jk) = ( .5*( zdiffus(ji-1,jj+1,jk) + zdiffus(ji-1,jj-1,jk)    &
1150                        &                   +zdiffus(ji+1,jj+1,jk) + zdiffus(ji+1,jj-1,jk) )  &
1151                        &              +1.*( zdiffus(ji-1,jj  ,jk) + zdiffus(ji  ,jj+1,jk)    &
1152                        &                   +zdiffus(ji  ,jj-1,jk) + zdiffus(ji+1,jj  ,jk) )  &
1153                        &              +2.*  zdiffus(ji  ,jj  ,jk)                          ) * etmean(ji,jj,jk) 
1154#endif               
1155                  END DO
1156               END DO
1157            END DO
1158         
1159         END SELECT
1160
1161         DO jk = 2, jpkm1                       ! vertical slab
1162            !
1163            !  Minimum value on the eddy diffusivity
1164            ! ----------------------------------------
1165            DO jj = 2, jpjm1
1166               DO ji = fs_2, fs_jpim1   ! vector opt.
1167                  avt(ji,jj,jk) = MAX( avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1168#if defined key_zdfddm 
1169                  avs(ji,jj,jk) = MAX( avs(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1170#endif
1171               END DO
1172            END DO
1173
1174            !
1175            ! Minimum value on the eddy viscosity
1176            ! ----------------------------------------
1177            DO jj = 1, jpj
1178               DO ji = 1, jpi
1179                  avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
1180                  avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
1181               END DO
1182            END DO
1183            !
1184         END DO
1185
1186         ! Lateral boundary conditions on avt  (sign unchanged)
1187         CALL lbc_lnk( hkpp(:,:), 'T', 1. )
1188
1189         ! Lateral boundary conditions on avt  (sign unchanged)
1190         CALL lbc_lnk( avt(:,:,:), 'W', 1. )
1191#if defined key_zdfddm 
1192         CALL lbc_lnk( avs(:,:,:), 'W', 1. ) 
1193#endif
1194         ! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
1195         CALL lbc_lnk( avmu(:,:,:), 'U', 1. )   ;    CALL lbc_lnk( avmv(:,:,:), 'V', 1. ) 
1196 
1197         IF(ln_ctl) THEN
1198#if defined key_zdfddm
1199            CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
1200#else
1201            CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', ovlap=1, kdim=jpk)
1202#endif
1203            CALL prt_ctl(tab3d_1=avmu, clinfo1=' kpp - u: ', mask1=umask,  &
1204               &         tab3d_2=avmv, clinfo2=      ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
1205         ENDIF
1206
1207      CALL wrk_dealloc( jpi, zmoa, zekman, zhmax, zria, zhbl )
1208      CALL wrk_dealloc( jpi, za2m, za3m, zkmpm, za2t, za3t, zkmpt )
1209      CALL wrk_dealloc( jpi,2, zriblk )
1210      CALL wrk_dealloc( jpi,3, zmoek, kjstart = 0 )
1211      CALL wrk_dealloc( jpi,3, zdept )
1212      CALL wrk_dealloc( jpi,4, zdepw, zdift, zvisc )
1213      CALL wrk_dealloc( jpi,jpj, zBo, zBosol, zustar )
1214      CALL wrk_dealloc( jpi,jpk, zmask, zblcm, zblct )
1215#if defined key_zdfddm
1216      CALL wrk_dealloc( jpi,4, zdifs )
1217      CALL wrk_dealloc( jpi, zmoa, za2s, za3s, zkmps )
1218      CALL wrk_dealloc( jpi,jpk, zblcs )
1219      CALL wrk_dealloc( jpi,jpi,jpk, zdiffus )
1220#endif
1221      !
1222      IF( nn_timing == 1 )  CALL timing_stop('zdf_kpp')
1223      !
1224   END SUBROUTINE zdf_kpp
1225
1226
1227   SUBROUTINE tra_kpp( kt )
1228      !!----------------------------------------------------------------------
1229      !!                  ***  ROUTINE tra_kpp  ***
1230      !!
1231      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
1232      !!
1233      !! ** Method  :   ???
1234      !!----------------------------------------------------------------------
1235      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
1236      !!----------------------------------------------------------------------
1237      INTEGER, INTENT(in) :: kt
1238      INTEGER :: ji, jj, jk
1239      !
1240      IF( nn_timing == 1 )  CALL timing_start('tra_kpp')
1241      !
1242      IF( kt == nit000 ) THEN
1243         IF(lwp) WRITE(numout,*) 
1244         IF(lwp) WRITE(numout,*) 'tra_kpp : KPP non-local tracer fluxes'
1245         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1246      ENDIF
1247
1248      IF( l_trdtra )   THEN                    !* Save ta and sa trends
1249         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
1250         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
1251      ENDIF
1252
1253      ! add non-local temperature and salinity flux ( in convective case only)
1254      DO jk = 1, jpkm1
1255         DO jj = 2, jpjm1 
1256            DO ji = fs_2, fs_jpim1
1257               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
1258                  &                 - (  ghats(ji,jj,jk  ) * avt  (ji,jj,jk  )   & 
1259                  &                    - ghats(ji,jj,jk+1) * avt  (ji,jj,jk+1) ) * wt0(ji,jj) / fse3t(ji,jj,jk)
1260               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
1261                  &                 - (  ghats(ji,jj,jk  ) * fsavs(ji,jj,jk  )   & 
1262                  &                    - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * ws0(ji,jj) / fse3t(ji,jj,jk)
1263            END DO
1264         END DO
1265      END DO
1266
1267      ! save the non-local tracer flux trends for diagnostic
1268      IF( l_trdtra )   THEN
1269         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
1270         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
1271!!bug gm jpttdzdf ==> jpttkpp
1272         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt )
1273         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds )
1274         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
1275      ENDIF
1276
1277      IF(ln_ctl) THEN
1278         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' kpp  - Ta: ', mask1=tmask,   &
1279         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
1280      ENDIF
1281      !
1282      IF( nn_timing == 1 )  CALL timing_stop('tra_kpp')
1283      !
1284   END SUBROUTINE tra_kpp
1285
1286#if defined key_top
1287   !!----------------------------------------------------------------------
1288   !!   'key_top'                                                TOP models
1289   !!----------------------------------------------------------------------
1290   SUBROUTINE trc_kpp( kt )
1291      !!----------------------------------------------------------------------
1292      !!                  ***  ROUTINE trc_kpp  ***
1293      !!
1294      !! ** Purpose :   compute and add to the tracer trend the non-local
1295      !!                tracer flux
1296      !!
1297      !! ** Method  :   ???
1298      !!
1299      !! history :
1300      !!            9.0  ! 2005-11 (G. Madec)  Original code
1301      !!       NEMO 3.3  ! 2010-06 (C. Ethe )  Adapted to passive tracers
1302      !!----------------------------------------------------------------------
1303      USE trc
1304      USE prtctl_trc          ! Print control
1305      !
1306      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
1307      !
1308      INTEGER  ::   ji, jj, jk, jn      ! Dummy loop indices
1309      CHARACTER (len=35) :: charout
1310      REAL(wp) ::   ztra, zflx
1311      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrtrd
1312      !!----------------------------------------------------------------------
1313
1314      IF( kt == nit000 ) THEN
1315         IF(lwp) WRITE(numout,*) 
1316         IF(lwp) WRITE(numout,*) 'trc_kpp : KPP non-local tracer fluxes'
1317         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1318      ENDIF
1319
1320      IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) )
1321      !
1322      DO jn = 1, jptra
1323         !
1324         IF( l_trdtrc )  ztrtrd(:,:,:)  = tra(:,:,:,jn)
1325         ! add non-local on passive tracer flux ( in convective case only)
1326         DO jk = 1, jpkm1
1327            DO jj = 2, jpjm1 
1328               DO ji = fs_2, fs_jpim1
1329                  ! Surface tracer flux for non-local term
1330                  zflx = - ( sfx (ji,jj) * tra(ji,jj,1,jn) * rcs ) * tmask(ji,jj,1)
1331                  ! compute the trend
1332                  ztra = - ( ghats(ji,jj,jk  ) * fsavs(ji,jj,jk  )   &
1333                  &        - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * zflx / fse3t(ji,jj,jk)
1334                  ! add the trend to the general trend
1335                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn)  + ztra
1336               END DO
1337            END DO
1338         END DO
1339         !
1340         IF( l_trdtrc ) THEN         ! save the non-local tracer flux trends for diagnostic
1341            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
1342            CALL trd_tra( kt, 'TRC', jn, jptra_trd_zdf, ztrtrd(:,:,:) )
1343         ENDIF
1344         !
1345      END DO
1346      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
1347      IF( ln_ctl )   THEN
1348         WRITE(charout, FMT="(' kpp')")  ;  CALL prt_ctl_trc_info(charout)
1349         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
1350      ENDIF
1351      !
1352   END SUBROUTINE trc_kpp
1353
1354#else
1355   !!----------------------------------------------------------------------
1356   !!   NO 'key_top'           DUMMY routine                  No TOP models
1357   !!----------------------------------------------------------------------
1358   SUBROUTINE trc_kpp( kt )         ! Dummy routine
1359      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
1360      WRITE(*,*) 'tra_kpp: You should not have seen this print! error?', kt
1361   END SUBROUTINE trc_kpp
1362#endif
1363
1364   SUBROUTINE zdf_kpp_init
1365      !!----------------------------------------------------------------------
1366      !!                  ***  ROUTINE zdf_kpp_init  ***
1367      !!                     
1368      !! ** Purpose :   Initialization of the vertical eddy diffivity and
1369      !!      viscosity when using a kpp turbulent closure scheme
1370      !!
1371      !! ** Method  :   Read the namkpp namelist and check the parameters
1372      !!      called at the first timestep (nit000)
1373      !!
1374      !! ** input   :   Namlist namkpp
1375      !!----------------------------------------------------------------------
1376      INTEGER  ::   ji, jj, jk     ! dummy loop indices
1377#if ! defined key_kppcustom
1378      INTEGER  ::   jm             ! dummy loop indices     
1379      REAL(wp) ::   zref, zdist    ! tempory scalars
1380#endif
1381#if defined key_kpplktb
1382      REAL(wp) ::   zustar, zucube, zustvk, zeta, zehat   ! tempory scalars
1383#endif
1384      INTEGER  ::   ios            ! Local integer output status for namelist read
1385      REAL(wp) ::   zhbf           ! tempory scalars
1386      LOGICAL  ::   ll_kppcustom   ! 1st ocean level taken as surface layer
1387      LOGICAL  ::   ll_kpplktb     ! Lookup table for turbul. velocity scales
1388      !!
1389      NAMELIST/namzdf_kpp/ ln_kpprimix, rn_difmiw, rn_difsiw, rn_riinfty, rn_difri, rn_bvsqcon, rn_difcon, nn_ave
1390      !!----------------------------------------------------------------------
1391      !
1392      IF( nn_timing == 1 )  CALL timing_start('zdf_kpp_init')
1393      !
1394      REWIND( numnam_ref )              ! Namelist namzdf_kpp in reference namelist : Vertical eddy diffivity and viscosity using kpp turbulent closure scheme
1395      READ  ( numnam_ref, namzdf_kpp, IOSTAT = ios, ERR = 901)
1396901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_kpp in reference namelist', lwp )
1397
1398      REWIND( numnam_cfg )              ! Namelist namzdf_kpp in configuration namelist : Vertical eddy diffivity and viscosity using kpp turbulent closure scheme
1399      READ  ( numnam_cfg, namzdf_kpp, IOSTAT = ios, ERR = 902 )
1400902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_kpp in configuration namelist', lwp )
1401      WRITE ( numond, namzdf_kpp )
1402
1403      IF(lwp) THEN                    ! Control print
1404         WRITE(numout,*)
1405         WRITE(numout,*) 'zdf_kpp_init : K-Profile Parameterisation'
1406         WRITE(numout,*) '~~~~~~~~~~~~'
1407         WRITE(numout,*) '   Namelist namzdf_kpp : set tke mixing parameters'
1408         WRITE(numout,*) '     Shear instability mixing                      ln_kpprimix = ', ln_kpprimix
1409         WRITE(numout,*) '     max. internal wave viscosity                  rn_difmiw   = ', rn_difmiw
1410         WRITE(numout,*) '     max. internal wave diffusivity                rn_difsiw   = ', rn_difsiw
1411         WRITE(numout,*) '     Richardson Number limit for shear instability rn_riinfty  = ', rn_riinfty
1412         WRITE(numout,*) '     max. shear mixing at Rig = 0                  rn_difri    = ', rn_difri
1413         WRITE(numout,*) '     Brunt-Vaisala squared for max. convection     rn_bvsqcon  = ', rn_bvsqcon
1414         WRITE(numout,*) '     max. mix. in interior convec.                 rn_difcon   = ', rn_difcon
1415         WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
1416      ENDIF
1417
1418      !                              ! allocate zdfkpp arrays
1419      IF( zdf_kpp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_kpp_init : unable to allocate arrays' )
1420
1421      ll_kppcustom = .FALSE.
1422      ll_kpplktb   = .FALSE.
1423
1424#if defined key_kppcustom
1425      ll_kppcustom = .TRUE.
1426#endif
1427#if defined key_kpplktb
1428      ll_kpplktb   = .TRUE.
1429#endif
1430      IF(lwp) THEN
1431         WRITE(numout,*) '     Lookup table for turbul. velocity scales ll_kpplktb   = ', ll_kpplktb
1432         WRITE(numout,*) '     1st ocean level taken as surface layer   ll_kppcustom = ', ll_kppcustom
1433      ENDIF
1434     
1435      IF( lk_zdfddm) THEN
1436         IF(lwp) THEN
1437            WRITE(numout,*)
1438            WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
1439            WRITE(numout,*) '    CAUTION : done in routine zdfkpp, not in routine zdfddm '
1440         ENDIF
1441      ENDIF
1442     
1443
1444      !set constants not in namelist
1445      !-----------------------------
1446      Vtc  = rconcv * SQRT( 0.2 / ( rconcs * epsilon ) ) / ( vonk * vonk * Ricr )
1447      rcg  = rcstar * vonk * ( rconcs * vonk * epsilon )**pthird
1448
1449      IF(lwp) THEN
1450         WRITE(numout,*)
1451         WRITE(numout,*) '     Constant value for unreso. turbul. velocity shear Vtc = ', Vtc
1452         WRITE(numout,*) '     Non-dimensional coef. for nonlocal transport      rcg = ', rcg
1453       ENDIF
1454
1455      ! ratt is the attenuation coefficient for solar flux
1456      ! Should be different is s_coordinate
1457      DO jk = 1, jpk
1458         zhbf     = - fsdept(1,1,jk) * hbf
1459         ratt(jk) = 1.0 - ( rabs * EXP( zhbf / xsi1 ) + ( 1.0 - rabs ) * EXP( zhbf / xsi2 ) )       
1460      ENDDO
1461
1462      ! Horizontal average : initialization of weighting arrays
1463      ! -------------------
1464     
1465      SELECT CASE ( nn_ave )
1466
1467      CASE ( 0 )                ! no horizontal average
1468         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
1469         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
1470         ! weighting mean arrays etmean, eumean and evmean
1471         !           ( 1  1 )                                          ( 1 )
1472         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
1473         !                         
1474         etmean(:,:,:) = 0.e0
1475         eumean(:,:,:) = 0.e0
1476         evmean(:,:,:) = 0.e0
1477         
1478         DO jk = 1, jpkm1
1479            DO jj = 2, jpjm1
1480               DO ji = 2, jpim1   ! vector opt.
1481                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
1482                  &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
1483                  &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
1484                 
1485                  eumean(ji,jj,jk) = umask(ji,jj,jk)                     &
1486                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk)  )
1487
1488                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                     &
1489                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk)  )
1490               END DO
1491            END DO
1492         END DO
1493
1494      CASE ( 1 )                ! horizontal average
1495         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
1496         ! weighting mean arrays etmean, eumean and evmean
1497         !           ( 1/2  1  1/2 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
1498         ! avt = 1/8 ( 1    2  1   )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
1499         !           ( 1/2  1  1/2 )              ( 1/2  1/2 )
1500         etmean(:,:,:) = 0.e0
1501         eumean(:,:,:) = 0.e0
1502         evmean(:,:,:) = 0.e0
1503         
1504         DO jk = 1, jpkm1
1505            DO jj = 2, jpjm1
1506               DO ji = fs_2, fs_jpim1   ! vector opt.
1507                  etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
1508                     & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
1509                     &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
1510                     &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
1511                     &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
1512                     &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
1513                 
1514                  eumean(ji,jj,jk) = umask(ji,jj,jk)                        &
1515                     &  / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
1516                     &       +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
1517                     &              +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
1518                 
1519                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                        &
1520                     &  / MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
1521                     &       +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
1522                     &              +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
1523               END DO
1524            END DO
1525         END DO
1526
1527      CASE DEFAULT
1528         WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
1529         CALL ctl_stop( ctmp1 )
1530
1531      END SELECT
1532 
1533      ! Initialization of vertical eddy coef. to the background value
1534      ! -------------------------------------------------------------
1535      DO jk = 1, jpk
1536         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1537         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1538         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1539      END DO
1540
1541      ! zero the surface flux for non local term and kpp mixed layer depth
1542      ! ------------------------------------------------------------------
1543      ghats(:,:,:) = 0.
1544      wt0  (:,:  ) = 0.
1545      ws0  (:,:  ) = 0.
1546      hkpp (:,:  ) = 0. ! just a diagnostic (not essential)
1547
1548#if ! defined key_kppcustom
1549      ! compute arrays (del, wz) for reference mean values
1550      ! (increase speed for vectorization key_kppcustom not defined)
1551      del(1:jpk, 1:jpk) = 0.
1552      DO jk = 1, jpk
1553         zref = epsilon * fsdept(1,1,jk)   
1554         DO jm = 1 , jpk
1555            zdist = zref - fsdepw(1,1,jm)   
1556            IF( zdist > 0.  ) THEN
1557               del(jk,jm) = MIN( zdist, fse3t(1,1,jm) ) / zref   
1558            ELSE
1559               del(jk,jm) = 0.
1560            ENDIF
1561         ENDDO
1562      ENDDO
1563#endif
1564
1565#if defined key_kpplktb
1566      ! build lookup table for turbulent velocity scales
1567      dezehat = ( dehatmax - dehatmin ) / nilktbm1
1568      deustar = ( ustmax   - ustmin   ) / njlktbm1
1569 
1570      DO jj = 1, njlktb
1571         zustar = ( jj - 1) * deustar + ustmin
1572         zustvk = vonk * zustar 
1573         zucube = zustar * zustar * zustar 
1574         DO ji = 1 , nilktb
1575            zehat = ( ji - 1 ) * dezehat + dehatmin
1576            zeta   = zehat / ( zucube + epsln )
1577            IF( zehat >= 0 ) THEN             ! Stable case
1578               wmlktb(ji,jj) = zustvk / ABS( 1.0 + rconc1 * zeta + epsln )                       
1579               wslktb(ji,jj) = wmlktb(ji,jj)
1580            ELSE                                ! Unstable case
1581               IF( zeta > rzetam ) THEN
1582                  wmlktb(ji,jj) = zustvk * ABS( 1.0    - rconc2 * zeta )**pfourth
1583               ELSE
1584                  wmlktb(ji,jj) = zustvk * ABS( rconam - rconcm * zeta )**pthird
1585               ENDIF
1586               
1587               IF( zeta > rzetas ) THEN
1588                  wslktb(ji,jj) = zustvk * SQRT( ABS( 1.0 - rconc2 * zeta ) )
1589               ELSE
1590                  wslktb(ji,jj) = zustvk * ABS( rconas - rconcs * zeta )**pthird
1591               ENDIF
1592            ENDIF
1593         END DO
1594      END DO
1595#endif
1596      !
1597      IF( nn_timing == 1 )  CALL timing_stop('zdf_kpp_init')
1598      !
1599   END SUBROUTINE zdf_kpp_init
1600
1601#else
1602   !!----------------------------------------------------------------------
1603   !!   Dummy module :                                        NO KPP scheme
1604   !!----------------------------------------------------------------------
1605   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfkpp = .FALSE.   !: KPP flag
1606CONTAINS
1607   SUBROUTINE zdf_kpp_init           ! Dummy routine
1608      WRITE(*,*) 'zdf_kpp_init: You should not have seen this print! error?'
1609   END SUBROUTINE zdf_kpp_init
1610   SUBROUTINE zdf_kpp( kt )          ! Dummy routine
1611      WRITE(*,*) 'zdf_kpp: You should not have seen this print! error?', kt
1612   END SUBROUTINE zdf_kpp
1613   SUBROUTINE tra_kpp( kt )          ! Dummy routine
1614      WRITE(*,*) 'tra_kpp: You should not have seen this print! error?', kt
1615   END SUBROUTINE tra_kpp
1616   SUBROUTINE trc_kpp( kt )          ! Dummy routine
1617      WRITE(*,*) 'trc_kpp: You should not have seen this print! error?', kt
1618   END SUBROUTINE trc_kpp
1619#endif
1620
1621   !!======================================================================
1622END MODULE zdfkpp
Note: See TracBrowser for help on using the repository browser.