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

Last change on this file since 3186 was 3186, checked in by smasson, 10 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo

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