source: trunk/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 888

Last change on this file since 888 was 888, checked in by ctlod, 13 years ago

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113

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