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

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

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 3764

Last change on this file since 3764 was 3764, checked in by smasson, 11 years ago

dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012

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