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

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

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 4990

Last change on this file since 4990 was 4990, checked in by timgraham, 9 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

  • Property svn:keywords set to Id
File size: 78.1 KB
RevLine 
[255]1MODULE zdfkpp
2   !!======================================================================
3   !!                       ***  MODULE  zdfkpp  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the KPP
5   !!                 turbulent closure parameterization
6   !!=====================================================================
[2528]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
[2715]11   !!            3.3  ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
[503]12   !!----------------------------------------------------------------------
[255]13#if defined key_zdfkpp   ||   defined key_esopa
14   !!----------------------------------------------------------------------
15   !!   'key_zdfkpp'                                             KPP scheme
16   !!----------------------------------------------------------------------
[3625]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)
[255]21   !!----------------------------------------------------------------------
[3625]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
[4990]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   !
[3625]33   USE in_out_manager ! I/O manager
34   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
35   USE prtctl         ! Print control
[4990]36   USE wrk_nemo       ! work arrays
[3625]37   USE timing         ! Timing
[4990]38   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[255]39
40   IMPLICIT NONE
41   PRIVATE
42
[2528]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
[255]47
[1537]48   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfkpp = .TRUE.    !: KPP vertical mixing flag
49
[2715]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
[503]54
[4147]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
[255]63
64#if defined key_zdfddm
[2528]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
[255]69#endif
[4147]70   LOGICAL  ::   ln_kpprimix       ! Shear instability mixing
[255]71
[2528]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
[255]76
[2528]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
[255]102
103#if ! defined key_kppcustom
[2715]104   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   del     ! array for reference mean values of vertical integration
[255]105#endif
106
107#if defined key_kpplktb
[2528]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   !
[255]113
[2528]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)
[255]116
[2528]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
[255]123#endif
[2715]124   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   ratt   ! attenuation coef  (already defines in module traqsr,
[1601]125   !                                    ! but only if the solar radiation penetration is considered)
[2528]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
[255]131      !                           ! (default values: water type Ib)
132
[2715]133   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean, eumean, evmean   ! coeff. used for hor. smoothing at t-, u- & v-points
[2528]134       
[899]135#if defined key_c1d
[2715]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
[255]141#endif
142
[1601]143   INTEGER  ::   jip = 62 , jjp = 111
[255]144
145   !! * Substitutions
146#  include "domzgr_substitute.h90"
147#  include "vectopt_loop_substitute.h90"
[896]148#  include  "zdfddm_substitute.h90"
[255]149   !!----------------------------------------------------------------------
[2715]150   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[888]151   !! $Id$
[2715]152   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[255]153   !!----------------------------------------------------------------------
154CONTAINS
155
[2715]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
[2528]177   SUBROUTINE zdf_kpp( kt )
[255]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      !!
[503]207      !! References : Large W.G., Mc Williams J.C. and Doney S.C.             
[255]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      !!----------------------------------------------------------------------
[2528]212      USE oce     , zviscos => ua   ! temp. array for viscosities use ua as workspace
[3294]213      USE oce     , zdiffut => va   ! temp. array for diffusivities use sa as workspace
[503]214      !!
215      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
216      !!
217      INTEGER ::   ji, jj, jk              ! dummy loop indices
218      INTEGER ::   ikbot, jkmax, jkm1, jkp2   !
[255]219
[2528]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
[255]225#if defined key_kpplktb
[2528]226      INTEGER ::    il, jl                                          ! Lookup table or Analytical functions
227      REAL(wp) ::   ud, zfrac, ufrac, zwam, zwbm, zwas, zwbs        !
[255]228#else
[2528]229      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm      !
[255]230#endif
[2528]231      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density
[255]232#if ! defined key_kppcustom     
[2528]233      INTEGER  ::   jm                          ! dummy loop indices
234      REAL(wp) ::   zr1, zr2, zr3, zr4, zrhop   ! Compression terms
[255]235#endif
[2528]236      REAL(wp) ::   zflag, ztemp, zrn2, zdep21, zdep32, zdep43
237      REAL(wp) ::   zdku2, zdkv2, ze3sqr, zsh2, zri, zfri          ! Interior richardson mixing
[3294]238      REAL(wp), POINTER, DIMENSION(:,:) ::   zmoek                 ! Moning-Obukov limitation
[2715]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 
[2528]245      REAL(wp) ::   zflagri, zflagek, zflagmo, zflagh, zflagkb   !
[2715]246      REAL(wp), POINTER, DIMENSION(:)   ::   za2m, za3m, zkmpm, za2t, za3t, zkmpt   ! Shape function (G)
[2528]247      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
[255]248#if defined key_zdfddm
[4990]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
[3764]252      REAL(wp), POINTER, DIMENSION(:)     ::   za2s, za3s, zkmps
253      REAL(wp), POINTER, DIMENSION(:,:)   ::   zblcs
[3294]254      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdiffus
[255]255#endif
[3294]256      REAL(wp), POINTER, DIMENSION(:,:) ::   zBo, zBosol, zustar         ! Surface buoyancy forcing, friction velocity
257      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask, zblcm, zblct
[255]258      !!--------------------------------------------------------------------
[3294]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 )
[3764]269      CALL wrk_alloc( jpi,jpk, zmask, zblcm, zblct )
[2715]270#if defined key_zdfddm
[3294]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 )
[2715]275#endif
276
[4990]277      zviscos(:,:,:) = 0._wp
278      zblcm  (:,:  ) = 0._wp
279      zdiffut(:,:,:) = 0._wp
280      zblct  (:,:  ) = 0._wp 
[255]281#if defined key_zdfddm
[4990]282      zdiffus(:,:,:) = 0._wp
283      zblcs  (:,:  ) = 0._wp 
[255]284#endif
[4990]285      ghats  (:,:,:) = 0._wp
286      zBo    (:,:  ) = 0._wp
287      zBosol (:,:  ) = 0._wp
288      zustar (:,:  ) = 0._wp
289      !
[255]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               ! -------------------------------------
[1537]298               avmu(ji,jj,jk)  = rn_difmiw 
299               avt (ji,jj,jk)  = rn_difsiw             
[255]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
[1082]317                  zsh2   = 0.5 * ( zdku2 + zdkv2 ) * ze3sqr
[255]318                  zri    = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + epsln ) 
[899]319#if defined key_c1d
[255]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. )
[1537]326                  zfri  = MIN( zfri / rn_riinfty , 1.0 )
[255]327                  zfri  = ( 1.0 - zfri * zfri )
328                  zfri  = zfri * zfri  * zfri
329                  ! add shear contribution to mixing coef.
[1537]330                  avmu(ji,jj,jk) =  avmu(ji,jj,jk) + rn_difri * zfri   
331                  avt (ji,jj,jk) =  avt (ji,jj,jk) + rn_difri * zfri   
[255]332               ENDIF
[4990]333               !
[255]334#if defined key_zdfddm 
[4990]335               !
[255]336               !  Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90
[4990]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
[255]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                  !
[4990]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);
[255]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
[4990]371                  IF( zrrau < 0.5) THEN   ;   zavdds = zavddt * 0.15 * zrrau
372                  ELSE                    ;   zavdds = zavddt * (1.85 * zrrau - 0.85 ) 
[255]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
[4990]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            !
[255]397            ! Radiative surface buoyancy force
398            zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
399            ! Non radiative surface buoyancy force
[3625]400            zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * ( emp(ji,jj)-rnf(ji,jj) )  &
[3792]401               &                                         -  grav * zbeta * rcs * sfx(ji,jj)
[255]402            ! Surface Temperature flux for non-local term
[3625]403            wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* r1_rau0_rcp * tmask(ji,jj,1)
[255]404            ! Surface salinity flux for non-local term
[3625]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) 
[4990]407         END DO
408      END DO
[255]409
[1601]410      zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. ) 
[255]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
[1695]418            zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos +  epsln ) )
[4990]419         END DO
420      END DO
[255]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
[2528]446            ikbot     = mbkt(ji,jj)    ! ikbot is the last T point in the water
[255]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 ) )   
[899]453#if defined key_c1d
[255]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)
[3294]487               zt   = tsn(ji,jj,1,jp_tem)
488               zs   = tsn(ji,jj,1,jp_sal)
[255]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
[3294]500                  zt   = zt  + del(jk,jm) * tsn(ji,jj,jm,jp_tem)
501                  zs   = zs  + del(jk,jm) * tsn(ji,jj,jm,jp_sal)
[255]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
[3294]511               zsr = SQRT( ABS( tsn(ji,jj,jk,jp_sal) ) )
[255]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
[1601]536               IF ( nn_eos < 1 ) THEN
[255]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 
[899]682#if defined key_c1d
[255]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
[3764]1054               zdiffus(ji,jj,jk) = ( 1.0 - zmask(ji,jk) )         * avs (ji,jj,jk) & ! interior diffusivities
[255]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
[1537]1074      SELECT CASE ( nn_ave )
[255]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 
[258]1167         IF(ln_ctl) THEN
1168#if defined key_zdfddm
[516]1169            CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
[258]1170#else
[516]1171            CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', ovlap=1, kdim=jpk)
[258]1172#endif
[516]1173            CALL prt_ctl(tab3d_1=avmu, clinfo1=' kpp - u: ', mask1=umask,  &
1174               &         tab3d_2=avmv, clinfo2=      ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
[258]1175         ENDIF
1176
[3294]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 )
[3764]1184      CALL wrk_dealloc( jpi,jpk, zmask, zblcm, zblct )
[3294]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
[2715]1191      !
[3294]1192      IF( nn_timing == 1 )  CALL timing_stop('zdf_kpp')
1193      !
[255]1194   END SUBROUTINE zdf_kpp
1195
1196
[896]1197   SUBROUTINE tra_kpp( kt )
[463]1198      !!----------------------------------------------------------------------
1199      !!                  ***  ROUTINE tra_kpp  ***
1200      !!
[2528]1201      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
[463]1202      !!
1203      !! ** Method  :   ???
1204      !!----------------------------------------------------------------------
[2528]1205      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
[463]1206      !!----------------------------------------------------------------------
[896]1207      INTEGER, INTENT(in) :: kt
1208      INTEGER :: ji, jj, jk
[3294]1209      !
1210      IF( nn_timing == 1 )  CALL timing_start('tra_kpp')
1211      !
[463]1212      IF( kt == nit000 ) THEN
[2528]1213         IF(lwp) WRITE(numout,*) 
[463]1214         IF(lwp) WRITE(numout,*) 'tra_kpp : KPP non-local tracer fluxes'
1215         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1216      ENDIF
1217
[2528]1218      IF( l_trdtra )   THEN                    !* Save ta and sa trends
1219         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
1220         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[463]1221      ENDIF
1222
1223      ! add non-local temperature and salinity flux ( in convective case only)
1224      DO jk = 1, jpkm1
[2528]1225         DO jj = 2, jpjm1 
[463]1226            DO ji = fs_2, fs_jpim1
[2528]1227               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
1228                  &                 - (  ghats(ji,jj,jk  ) * avt  (ji,jj,jk  )   & 
1229                  &                    - ghats(ji,jj,jk+1) * avt  (ji,jj,jk+1) ) * wt0(ji,jj) / fse3t(ji,jj,jk)
1230               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
1231                  &                 - (  ghats(ji,jj,jk  ) * fsavs(ji,jj,jk  )   & 
1232                  &                    - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * ws0(ji,jj) / fse3t(ji,jj,jk)
[463]1233            END DO
1234         END DO
1235      END DO
1236
1237      ! save the non-local tracer flux trends for diagnostic
1238      IF( l_trdtra )   THEN
[2528]1239         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
1240         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
[463]1241!!bug gm jpttdzdf ==> jpttkpp
[4990]1242         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
1243         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
[2528]1244         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
[463]1245      ENDIF
1246
[2528]1247      IF(ln_ctl) THEN
1248         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' kpp  - Ta: ', mask1=tmask,   &
1249         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
[463]1250      ENDIF
[3294]1251      !
1252      IF( nn_timing == 1 )  CALL timing_stop('tra_kpp')
1253      !
[463]1254   END SUBROUTINE tra_kpp
1255
[2528]1256#if defined key_top
1257   !!----------------------------------------------------------------------
1258   !!   'key_top'                                                TOP models
1259   !!----------------------------------------------------------------------
1260   SUBROUTINE trc_kpp( kt )
1261      !!----------------------------------------------------------------------
1262      !!                  ***  ROUTINE trc_kpp  ***
1263      !!
1264      !! ** Purpose :   compute and add to the tracer trend the non-local
1265      !!                tracer flux
1266      !!
1267      !! ** Method  :   ???
1268      !!
1269      !! history :
1270      !!            9.0  ! 2005-11 (G. Madec)  Original code
1271      !!       NEMO 3.3  ! 2010-06 (C. Ethe )  Adapted to passive tracers
1272      !!----------------------------------------------------------------------
1273      USE trc
1274      USE prtctl_trc          ! Print control
[2715]1275      !
1276      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
1277      !
[2528]1278      INTEGER  ::   ji, jj, jk, jn      ! Dummy loop indices
[3792]1279      CHARACTER (len=35) :: charout
[2528]1280      REAL(wp) ::   ztra, zflx
1281      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrtrd
1282      !!----------------------------------------------------------------------
[463]1283
[2528]1284      IF( kt == nit000 ) THEN
1285         IF(lwp) WRITE(numout,*) 
1286         IF(lwp) WRITE(numout,*) 'trc_kpp : KPP non-local tracer fluxes'
1287         IF(lwp) WRITE(numout,*) '~~~~~~~   '
1288      ENDIF
1289
1290      IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) )
1291      !
1292      DO jn = 1, jptra
1293         !
1294         IF( l_trdtrc )  ztrtrd(:,:,:)  = tra(:,:,:,jn)
1295         ! add non-local on passive tracer flux ( in convective case only)
1296         DO jk = 1, jpkm1
1297            DO jj = 2, jpjm1 
1298               DO ji = fs_2, fs_jpim1
1299                  ! Surface tracer flux for non-local term
[3625]1300                  zflx = - ( sfx (ji,jj) * tra(ji,jj,1,jn) * rcs ) * tmask(ji,jj,1)
[2528]1301                  ! compute the trend
1302                  ztra = - ( ghats(ji,jj,jk  ) * fsavs(ji,jj,jk  )   &
1303                  &        - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * zflx / fse3t(ji,jj,jk)
1304                  ! add the trend to the general trend
1305                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn)  + ztra
1306               END DO
1307            END DO
1308         END DO
1309         !
[3792]1310         IF( l_trdtrc ) THEN         ! save the non-local tracer flux trends for diagnostic
1311            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
[4990]1312            CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:) )
[3792]1313         ENDIF
1314         !
[2528]1315      END DO
1316      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
1317      IF( ln_ctl )   THEN
1318         WRITE(charout, FMT="(' kpp')")  ;  CALL prt_ctl_trc_info(charout)
[3792]1319         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
[2528]1320      ENDIF
1321      !
1322   END SUBROUTINE trc_kpp
[2715]1323
1324#else
1325   !!----------------------------------------------------------------------
1326   !!   NO 'key_top'           DUMMY routine                  No TOP models
1327   !!----------------------------------------------------------------------
1328   SUBROUTINE trc_kpp( kt )         ! Dummy routine
1329      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
1330      WRITE(*,*) 'tra_kpp: You should not have seen this print! error?', kt
1331   END SUBROUTINE trc_kpp
[2528]1332#endif
1333
[255]1334   SUBROUTINE zdf_kpp_init
1335      !!----------------------------------------------------------------------
1336      !!                  ***  ROUTINE zdf_kpp_init  ***
1337      !!                     
1338      !! ** Purpose :   Initialization of the vertical eddy diffivity and
1339      !!      viscosity when using a kpp turbulent closure scheme
1340      !!
1341      !! ** Method  :   Read the namkpp namelist and check the parameters
1342      !!      called at the first timestep (nit000)
1343      !!
1344      !! ** input   :   Namlist namkpp
1345      !!----------------------------------------------------------------------
[2528]1346      INTEGER  ::   ji, jj, jk     ! dummy loop indices
[4990]1347      INTEGER  ::   ios            ! local integer
[255]1348#if ! defined key_kppcustom
[2528]1349      INTEGER  ::   jm             ! dummy loop indices     
1350      REAL(wp) ::   zref, zdist    ! tempory scalars
[255]1351#endif
1352#if defined key_kpplktb
[2528]1353      REAL(wp) ::   zustar, zucube, zustvk, zeta, zehat   ! tempory scalars
[255]1354#endif
[2528]1355      REAL(wp) ::   zhbf           ! tempory scalars
1356      LOGICAL  ::   ll_kppcustom   ! 1st ocean level taken as surface layer
1357      LOGICAL  ::   ll_kpplktb     ! Lookup table for turbul. velocity scales
[1537]1358      !!
[1601]1359      NAMELIST/namzdf_kpp/ ln_kpprimix, rn_difmiw, rn_difsiw, rn_riinfty, rn_difri, rn_bvsqcon, rn_difcon, nn_ave
[255]1360      !!----------------------------------------------------------------------
[3294]1361      !
1362      IF( nn_timing == 1 )  CALL timing_start('zdf_kpp_init')
1363      !
[4147]1364      REWIND( numnam_ref )              ! Namelist namzdf_kpp in reference namelist : Vertical eddy diffivity and viscosity using kpp turbulent closure scheme
1365      READ  ( numnam_ref, namzdf_kpp, IOSTAT = ios, ERR = 901)
1366901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_kpp in reference namelist', lwp )
[255]1367
[4147]1368      REWIND( numnam_cfg )              ! Namelist namzdf_kpp in configuration namelist : Vertical eddy diffivity and viscosity using kpp turbulent closure scheme
1369      READ  ( numnam_cfg, namzdf_kpp, IOSTAT = ios, ERR = 902 )
1370902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_kpp in configuration namelist', lwp )
[4624]1371      IF(lwm) WRITE ( numond, namzdf_kpp )
[4147]1372
[1537]1373      IF(lwp) THEN                    ! Control print
[255]1374         WRITE(numout,*)
[1537]1375         WRITE(numout,*) 'zdf_kpp_init : K-Profile Parameterisation'
[255]1376         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]1377         WRITE(numout,*) '   Namelist namzdf_kpp : set tke mixing parameters'
[1537]1378         WRITE(numout,*) '     Shear instability mixing                      ln_kpprimix = ', ln_kpprimix
1379         WRITE(numout,*) '     max. internal wave viscosity                  rn_difmiw   = ', rn_difmiw
1380         WRITE(numout,*) '     max. internal wave diffusivity                rn_difsiw   = ', rn_difsiw
1381         WRITE(numout,*) '     Richardson Number limit for shear instability rn_riinfty  = ', rn_riinfty
1382         WRITE(numout,*) '     max. shear mixing at Rig = 0                  rn_difri    = ', rn_difri
1383         WRITE(numout,*) '     Brunt-Vaisala squared for max. convection     rn_bvsqcon  = ', rn_bvsqcon
1384         WRITE(numout,*) '     max. mix. in interior convec.                 rn_difcon   = ', rn_difcon
1385         WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
[255]1386      ENDIF
1387
[2715]1388      !                              ! allocate zdfkpp arrays
1389      IF( zdf_kpp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_kpp_init : unable to allocate arrays' )
1390
[255]1391      ll_kppcustom = .FALSE.
1392      ll_kpplktb   = .FALSE.
1393
1394#if defined key_kppcustom
1395      ll_kppcustom = .TRUE.
1396#endif
1397#if defined key_kpplktb
1398      ll_kpplktb   = .TRUE.
1399#endif
1400      IF(lwp) THEN
1401         WRITE(numout,*) '     Lookup table for turbul. velocity scales ll_kpplktb   = ', ll_kpplktb
1402         WRITE(numout,*) '     1st ocean level taken as surface layer   ll_kppcustom = ', ll_kppcustom
1403      ENDIF
1404     
1405      IF( lk_zdfddm) THEN
1406         IF(lwp) THEN
1407            WRITE(numout,*)
1408            WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
1409            WRITE(numout,*) '    CAUTION : done in routine zdfkpp, not in routine zdfddm '
1410         ENDIF
1411      ENDIF
1412     
1413
1414      !set constants not in namelist
1415      !-----------------------------
1416      Vtc  = rconcv * SQRT( 0.2 / ( rconcs * epsilon ) ) / ( vonk * vonk * Ricr )
1417      rcg  = rcstar * vonk * ( rconcs * vonk * epsilon )**pthird
1418
1419      IF(lwp) THEN
[1537]1420         WRITE(numout,*)
[255]1421         WRITE(numout,*) '     Constant value for unreso. turbul. velocity shear Vtc = ', Vtc
1422         WRITE(numout,*) '     Non-dimensional coef. for nonlocal transport      rcg = ', rcg
1423       ENDIF
1424
1425      ! ratt is the attenuation coefficient for solar flux
1426      ! Should be different is s_coordinate
1427      DO jk = 1, jpk
1428         zhbf     = - fsdept(1,1,jk) * hbf
1429         ratt(jk) = 1.0 - ( rabs * EXP( zhbf / xsi1 ) + ( 1.0 - rabs ) * EXP( zhbf / xsi2 ) )       
1430      ENDDO
1431
1432      ! Horizontal average : initialization of weighting arrays
1433      ! -------------------
1434     
[1537]1435      SELECT CASE ( nn_ave )
[255]1436
1437      CASE ( 0 )                ! no horizontal average
1438         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv'
1439         IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
1440         ! weighting mean arrays etmean, eumean and evmean
1441         !           ( 1  1 )                                          ( 1 )
1442         ! avt = 1/4 ( 1  1 )     avmu = 1/2 ( 1  1 )       avmv=  1/2 ( 1 )
1443         !                         
1444         etmean(:,:,:) = 0.e0
1445         eumean(:,:,:) = 0.e0
1446         evmean(:,:,:) = 0.e0
1447         
1448         DO jk = 1, jpkm1
1449            DO jj = 2, jpjm1
1450               DO ji = 2, jpim1   ! vector opt.
1451                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
1452                  &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
1453                  &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
1454                 
1455                  eumean(ji,jj,jk) = umask(ji,jj,jk)                     &
1456                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk)  )
1457
1458                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                     &
1459                  &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk)  )
1460               END DO
1461            END DO
1462         END DO
1463
1464      CASE ( 1 )                ! horizontal average
1465         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv'
1466         ! weighting mean arrays etmean, eumean and evmean
1467         !           ( 1/2  1  1/2 )              ( 1/2  1/2 )             ( 1/2  1  1/2 )
1468         ! avt = 1/8 ( 1    2  1   )   avmu = 1/4 ( 1    1   )   avmv= 1/4 ( 1/2  1  1/2 )
1469         !           ( 1/2  1  1/2 )              ( 1/2  1/2 )
1470         etmean(:,:,:) = 0.e0
1471         eumean(:,:,:) = 0.e0
1472         evmean(:,:,:) = 0.e0
1473         
1474         DO jk = 1, jpkm1
1475            DO jj = 2, jpjm1
1476               DO ji = fs_2, fs_jpim1   ! vector opt.
1477                  etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
1478                     & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
1479                     &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
1480                     &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
1481                     &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
1482                     &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
1483                 
1484                  eumean(ji,jj,jk) = umask(ji,jj,jk)                        &
1485                     &  / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   &
1486                     &       +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   &
1487                     &              +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  )
1488                 
1489                  evmean(ji,jj,jk) = vmask(ji,jj,jk)                        &
1490                     &  / MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
1491                     &       +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   &
1492                     &              +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  )
1493               END DO
1494            END DO
1495         END DO
1496
1497      CASE DEFAULT
[1537]1498         WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
[896]1499         CALL ctl_stop( ctmp1 )
[255]1500
1501      END SELECT
1502 
1503      ! Initialization of vertical eddy coef. to the background value
1504      ! -------------------------------------------------------------
1505      DO jk = 1, jpk
1506         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1507         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1508         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1509      END DO
1510
1511      ! zero the surface flux for non local term and kpp mixed layer depth
1512      ! ------------------------------------------------------------------
1513      ghats(:,:,:) = 0.
1514      wt0  (:,:  ) = 0.
1515      ws0  (:,:  ) = 0.
1516      hkpp (:,:  ) = 0. ! just a diagnostic (not essential)
1517
1518#if ! defined key_kppcustom
1519      ! compute arrays (del, wz) for reference mean values
1520      ! (increase speed for vectorization key_kppcustom not defined)
1521      del(1:jpk, 1:jpk) = 0.
1522      DO jk = 1, jpk
1523         zref = epsilon * fsdept(1,1,jk)   
1524         DO jm = 1 , jpk
1525            zdist = zref - fsdepw(1,1,jm)   
1526            IF( zdist > 0.  ) THEN
1527               del(jk,jm) = MIN( zdist, fse3t(1,1,jm) ) / zref   
1528            ELSE
1529               del(jk,jm) = 0.
1530            ENDIF
1531         ENDDO
1532      ENDDO
1533#endif
1534
1535#if defined key_kpplktb
1536      ! build lookup table for turbulent velocity scales
1537      dezehat = ( dehatmax - dehatmin ) / nilktbm1
1538      deustar = ( ustmax   - ustmin   ) / njlktbm1
1539 
1540      DO jj = 1, njlktb
1541         zustar = ( jj - 1) * deustar + ustmin
1542         zustvk = vonk * zustar 
1543         zucube = zustar * zustar * zustar 
1544         DO ji = 1 , nilktb
1545            zehat = ( ji - 1 ) * dezehat + dehatmin
1546            zeta   = zehat / ( zucube + epsln )
1547            IF( zehat >= 0 ) THEN             ! Stable case
1548               wmlktb(ji,jj) = zustvk / ABS( 1.0 + rconc1 * zeta + epsln )                       
1549               wslktb(ji,jj) = wmlktb(ji,jj)
1550            ELSE                                ! Unstable case
1551               IF( zeta > rzetam ) THEN
1552                  wmlktb(ji,jj) = zustvk * ABS( 1.0    - rconc2 * zeta )**pfourth
1553               ELSE
1554                  wmlktb(ji,jj) = zustvk * ABS( rconam - rconcm * zeta )**pthird
1555               ENDIF
1556               
1557               IF( zeta > rzetas ) THEN
1558                  wslktb(ji,jj) = zustvk * SQRT( ABS( 1.0 - rconc2 * zeta ) )
1559               ELSE
1560                  wslktb(ji,jj) = zustvk * ABS( rconas - rconcs * zeta )**pthird
1561               ENDIF
1562            ENDIF
1563         END DO
1564      END DO
1565#endif
[2715]1566      !
[3294]1567      IF( nn_timing == 1 )  CALL timing_stop('zdf_kpp_init')
1568      !
[255]1569   END SUBROUTINE zdf_kpp_init
1570
1571#else
1572   !!----------------------------------------------------------------------
1573   !!   Dummy module :                                        NO KPP scheme
1574   !!----------------------------------------------------------------------
1575   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfkpp = .FALSE.   !: KPP flag
1576CONTAINS
[2528]1577   SUBROUTINE zdf_kpp_init           ! Dummy routine
1578      WRITE(*,*) 'zdf_kpp_init: You should not have seen this print! error?'
1579   END SUBROUTINE zdf_kpp_init
1580   SUBROUTINE zdf_kpp( kt )          ! Dummy routine
[255]1581      WRITE(*,*) 'zdf_kpp: You should not have seen this print! error?', kt
1582   END SUBROUTINE zdf_kpp
[2528]1583   SUBROUTINE tra_kpp( kt )          ! Dummy routine
[463]1584      WRITE(*,*) 'tra_kpp: You should not have seen this print! error?', kt
1585   END SUBROUTINE tra_kpp
[2528]1586   SUBROUTINE trc_kpp( kt )          ! Dummy routine
1587      WRITE(*,*) 'trc_kpp: You should not have seen this print! error?', kt
1588   END SUBROUTINE trc_kpp
[255]1589#endif
1590
1591   !!======================================================================
1592END MODULE zdfkpp
Note: See TracBrowser for help on using the repository browser.