source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 4990

Last change on this file since 4990 was 4990, checked in by timgraham, 6 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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