source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 3229

Last change on this file since 3229 was 3229, checked in by charris, 10 years ago

Added timing calls to most significant routines in LDF, SBC and ZDF.

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