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

Last change on this file since 896 was 896, checked in by rblod, 13 years ago

Update KPP according to step reorganization (nemo_v2) …, see ticket #101

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