source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

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