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.
zdftke.F90 in branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

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

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

File size: 51.8 KB
RevLine 
[1531]1MODULE zdftke
[1239]2   !!======================================================================
[1531]3   !!                       ***  MODULE  zdftke  ***
[1239]4   !! Ocean physics:  vertical mixing coefficient computed from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
[1492]7   !! History :  OPA  !  1991-03  (b. blanke)  Original code
8   !!            7.0  !  1991-11  (G. Madec)   bug fix
9   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
10   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
11   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
12   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
13   !!            8.0  !  1997-07  (G. Madec)   lbc
14   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
15   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
16   !!             -   !  2004-10  (C. Ethe )  1D configuration
17   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
18   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
19   !!                 !           - tke penetration (wind steering)
20   !!                 !           - suface condition for tke & mixing length
21   !!                 !           - Langmuir cells
22   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
23   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
24   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
25   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
26   !!                 !                                + cleaning of the parameters + bugs correction
[2528]27   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
[5120]28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
[1239]29   !!----------------------------------------------------------------------
[1531]30#if defined key_zdftke   ||   defined key_esopa
[1239]31   !!----------------------------------------------------------------------
[1531]32   !!   'key_zdftke'                                   TKE vertical physics
[1239]33   !!----------------------------------------------------------------------
[3625]34   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
35   !!   tke_tke       : tke time stepping: update tke at now time step (en)
36   !!   tke_avn       : compute mixing length scale and deduce avm and avt
37   !!   zdf_tke_init  : initialization, namelist read, and parameters control
38   !!   tke_rst       : read/write tke restart in ocean restart file
[1239]39   !!----------------------------------------------------------------------
[2528]40   USE oce            ! ocean: dynamics and active tracers variables
41   USE phycst         ! physical constants
42   USE dom_oce        ! domain: ocean
43   USE domvvl         ! domain: variable volume layer
[1492]44   USE sbc_oce        ! surface boundary condition: ocean
[2528]45   USE zdf_oce        ! vertical physics: ocean variables
46   USE zdfmxl         ! vertical physics: mixed layer
[1492]47   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
48   USE prtctl         ! Print control
49   USE in_out_manager ! I/O manager
50   USE iom            ! I/O manager library
[2715]51   USE lib_mpp        ! MPP library
[3294]52   USE wrk_nemo       ! work arrays
53   USE timing         ! Timing
[3625]54   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[6487]55#if defined key_agrif
56   USE agrif_opa_interp
57   USE agrif_opa_update
58#endif
[1239]59
[6487]60
61
[1239]62   IMPLICIT NONE
63   PRIVATE
64
[2528]65   PUBLIC   zdf_tke        ! routine called in step module
66   PUBLIC   zdf_tke_init   ! routine called in opa module
67   PUBLIC   tke_rst        ! routine called in step module
[1239]68
[2715]69   LOGICAL , PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
[1239]70
[4147]71   !                      !!** Namelist  namzdf_tke  **
72   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
73   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
74   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
75   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
76   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
77   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
78   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
79   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
80   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
81   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
82   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
83   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
84   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean
[6491]85   REAL(wp) ::   rn_c      ! fraction of TKE added within the mixed layer by nn_etau
[4147]86   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
87   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
[1239]88
[4147]89   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
90   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
[6491]91   REAL(wp) ::   rhtau                     ! coefficient to relate MLD to htau when nn_htau == 2
[2528]92   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
93   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
[1239]94
[2715]95   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
[6491]96   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
97   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4
[2715]98   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
99#if defined key_c1d
100   !                                                                        !!** 1D cfg only  **   ('key_c1d')
101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales
102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
103#endif
[1492]104
[1239]105   !! * Substitutions
106#  include "domzgr_substitute.h90"
107#  include "vectopt_loop_substitute.h90"
108   !!----------------------------------------------------------------------
[2715]109   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2528]110   !! $Id$
111   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[1239]112   !!----------------------------------------------------------------------
113CONTAINS
114
[2715]115   INTEGER FUNCTION zdf_tke_alloc()
116      !!----------------------------------------------------------------------
117      !!                ***  FUNCTION zdf_tke_alloc  ***
118      !!----------------------------------------------------------------------
119      ALLOCATE(                                                                    &
[6491]120         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &     
[2715]121#if defined key_c1d
122         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
123         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
124#endif
[6487]125         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc      )
[2715]126         !
127      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
128      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
129      !
130   END FUNCTION zdf_tke_alloc
131
132
[1531]133   SUBROUTINE zdf_tke( kt )
[1239]134      !!----------------------------------------------------------------------
[1531]135      !!                   ***  ROUTINE zdf_tke  ***
[1239]136      !!
137      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
[1492]138      !!              coefficients using a turbulent closure scheme (TKE).
[1239]139      !!
[1492]140      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
141      !!              is computed from a prognostic equation :
142      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
143      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
144      !!                  + avt N^2                      ! stratif. destruc.
145      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
[1239]146      !!      with the boundary conditions:
[1695]147      !!         surface: en = max( rn_emin0, rn_ebb * taum )
[1239]148      !!         bottom : en = rn_emin
[1492]149      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
150      !!
151      !!        The now Turbulent kinetic energy is computed using the following
152      !!      time stepping: implicit for vertical diffusion term, linearized semi
153      !!      implicit for kolmogoroff dissipation term, and explicit forward for
154      !!      both buoyancy and shear production terms. Therefore a tridiagonal
155      !!      linear system is solved. Note that buoyancy and shear terms are
156      !!      discretized in a energy conserving form (Bruchard 2002).
157      !!
158      !!        The dissipative and mixing length scale are computed from en and
159      !!      the stratification (see tke_avn)
160      !!
161      !!        The now vertical eddy vicosity and diffusivity coefficients are
162      !!      given by:
163      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
164      !!              avt = max( avmb, pdl * avm                 ) 
[1239]165      !!              eav = max( avmb, avm )
[1492]166      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
167      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
[1239]168      !!
169      !! ** Action  :   compute en (now turbulent kinetic energy)
170      !!                update avt, avmu, avmv (before vertical eddy coef.)
171      !!
172      !! References : Gaspar et al., JGR, 1990,
173      !!              Blanke and Delecluse, JPO, 1991
174      !!              Mellor and Blumberg, JPO 2004
175      !!              Axell, JGR, 2002
[1492]176      !!              Bruchard OM 2002
[1239]177      !!----------------------------------------------------------------------
[1492]178      INTEGER, INTENT(in) ::   kt   ! ocean time step
179      !!----------------------------------------------------------------------
[1481]180      !
[3632]181      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
182         avt (:,:,:) = avt_k (:,:,:) 
183         avm (:,:,:) = avm_k (:,:,:) 
184         avmu(:,:,:) = avmu_k(:,:,:) 
185         avmv(:,:,:) = avmv_k(:,:,:) 
186      ENDIF 
187      !
[2528]188      CALL tke_tke      ! now tke (en)
[1492]189      !
[2528]190      CALL tke_avn      ! now avt, avm, avmu, avmv
191      !
[3632]192      avt_k (:,:,:) = avt (:,:,:) 
193      avm_k (:,:,:) = avm (:,:,:) 
194      avmu_k(:,:,:) = avmu(:,:,:) 
195      avmv_k(:,:,:) = avmv(:,:,:) 
196      !
[6487]197#if defined key_agrif
198      ! Update child grid f => parent grid
199      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
200#endif     
201     !
[1531]202   END SUBROUTINE zdf_tke
[1239]203
[1492]204
[1481]205   SUBROUTINE tke_tke
[1239]206      !!----------------------------------------------------------------------
[1492]207      !!                   ***  ROUTINE tke_tke  ***
208      !!
209      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
210      !!
211      !! ** Method  : - TKE surface boundary condition
[2528]212      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
[1492]213      !!              - source term due to shear (saved in avmu, avmv arrays)
214      !!              - Now TKE : resolution of the TKE equation by inverting
215      !!                a tridiagonal linear system by a "methode de chasse"
216      !!              - increase TKE due to surface and internal wave breaking
217      !!
218      !! ** Action  : - en : now turbulent kinetic energy)
219      !!              - avmu, avmv : production of TKE by shear at u and v-points
220      !!                (= Kz dz[Ub] * dz[Un] )
[1239]221      !! ---------------------------------------------------------------------
[1705]222      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
[2528]223!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
224!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
[1705]225      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
226      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
227      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
228      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
229      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
230      REAL(wp) ::   ztau  , zdif                    !    -         -
231      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
232      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
[2528]233!!bfr      REAL(wp) ::   zebot                           !    -         -
[3294]234      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
235      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc
236      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
[1239]237      !!--------------------------------------------------------------------
[1492]238      !
[3294]239      IF( nn_timing == 1 )  CALL timing_start('tke_tke')
240      !
241      CALL wrk_alloc( jpi,jpj, imlc )    ! integer
242      CALL wrk_alloc( jpi,jpj, zhlc ) 
243      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
244      !
[1695]245      zbbrau = rn_ebb / rau0       ! Local constant initialisation
[2528]246      zfact1 = -.5_wp * rdt 
247      zfact2 = 1.5_wp * rdt * rn_ediss
248      zfact3 = 0.5_wp       * rn_ediss
[1492]249      !
[5120]250      !
[1492]251      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
252      !                     !  Surface boundary condition on tke
253      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[5120]254      IF ( ln_isfcav ) THEN
255         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
256            DO ji = fs_2, fs_jpim1   ! vector opt.
257               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1)
258            END DO
259         END DO
260      END IF
[1695]261      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
[1481]262         DO ji = fs_2, fs_jpim1   ! vector opt.
[5120]263            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
[1481]264         END DO
265      END DO
[2528]266     
267!!bfr   - start commented area
[1492]268      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
269      !                     !  Bottom boundary condition on tke
270      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1719]271      !
272      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
273      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
274      ! The condition is coded here for completion but commented out until there is proof that the
275      ! computational cost is justified
276      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
277      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
[1662]278!CDIR NOVERRCHK
[1719]279!!    DO jj = 2, jpjm1
[1662]280!CDIR NOVERRCHK
[1719]281!!       DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]282!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
283!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
284!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
285!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
[1719]286!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
[2528]287!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
[1719]288!!       END DO
289!!    END DO
[2528]290!!bfr   - end commented area
[1492]291      !
292      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[2528]293      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
[1492]294         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[1239]295         !
[1492]296         !                        !* total energy produce by LC : cumulative sum over jk
[2528]297         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
[1239]298         DO jk = 2, jpk
[2528]299            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
[1239]300         END DO
[1492]301         !                        !* finite Langmuir Circulation depth
[1705]302         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
[2528]303         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
[1239]304         DO jk = jpkm1, 2, -1
[1492]305            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
306               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
[1705]307                  zus  = zcof * taum(ji,jj)
[1239]308                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
309               END DO
310            END DO
311         END DO
[1492]312         !                               ! finite LC depth
313         DO jj = 1, jpj 
[1239]314            DO ji = 1, jpi
315               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
316            END DO
317         END DO
[1705]318         zcof = 0.016 / SQRT( zrhoa * zcdrag )
[5120]319!CDIR NOVERRCHK
[1492]320         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
[5120]321!CDIR NOVERRCHK
[1239]322            DO jj = 2, jpjm1
[5120]323!CDIR NOVERRCHK
[1239]324               DO ji = fs_2, fs_jpim1   ! vector opt.
[1705]325                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
[1492]326                  !                                           ! vertical velocity due to LC
[1239]327                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
328                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
[1492]329                  !                                           ! TKE Langmuir circulation source term
[6487]330                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   &
331                     &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1239]332               END DO
333            END DO
334         END DO
335         !
336      ENDIF
[1492]337      !
338      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
339      !                     !  Now Turbulent kinetic energy (output in en)
340      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
341      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
342      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
343      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
344      !
345      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
346         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
347            DO ji = 1, jpi
348               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
349                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
[5120]350                  &                            / (  fse3uw_n(ji,jj,jk)               &
351                  &                              *  fse3uw_b(ji,jj,jk)  )
[1492]352               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
353                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
354                  &                            / (  fse3vw_n(ji,jj,jk)               &
355                  &                              *  fse3vw_b(ji,jj,jk)  )
356            END DO
357         END DO
358      END DO
359      !
[5120]360      DO jk = 2, jpkm1           !* Matrix and right hand side in en
361         DO jj = 2, jpjm1
362            DO ji = fs_2, fs_jpim1   ! vector opt.
[1492]363               zcof   = zfact1 * tmask(ji,jj,jk)
[6498]364# if defined key_zdftmx_new
365               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
366               zzd_up = zcof * ( MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) )   &  ! upper diagonal
367                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
368               zzd_lw = zcof * ( MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) )   &  ! lower diagonal
369                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
370# else
[1492]371               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
372                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
373               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
374                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
[6498]375# endif
[1492]376                  !                                                           ! shear prod. at w-point weightened by mask
[2528]377               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
378                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
[1492]379                  !
380               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
381               zd_lw(ji,jj,jk) = zzd_lw
[2528]382               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
[1239]383               !
[1492]384               !                                   ! right hand side in en
[1481]385               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
[4990]386                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
[5120]387                  &                                 * wmask(ji,jj,jk)
[1239]388            END DO
[5120]389         END DO
390      END DO
391      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
392      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
393         DO jj = 2, jpjm1
394            DO ji = fs_2, fs_jpim1    ! vector opt.
[1492]395               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
[1239]396            END DO
[5120]397         END DO
398      END DO
399      !
400      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
401      DO jj = 2, jpjm1
402         DO ji = fs_2, fs_jpim1   ! vector opt.
403            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
404         END DO
405      END DO
406      DO jk = 3, jpkm1
407         DO jj = 2, jpjm1
408            DO ji = fs_2, fs_jpim1    ! vector opt.
[1492]409               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
[1239]410            END DO
[5120]411         END DO
412      END DO
413      !
414      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
415      DO jj = 2, jpjm1
416         DO ji = fs_2, fs_jpim1   ! vector opt.
[1492]417            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
[5120]418         END DO
419      END DO
420      DO jk = jpk-2, 2, -1
421         DO jj = 2, jpjm1
422            DO ji = fs_2, fs_jpim1    ! vector opt.
[1492]423               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
[1239]424            END DO
[5120]425         END DO
426      END DO
427      DO jk = 2, jpkm1                             ! set the minimum value of tke
428         DO jj = 2, jpjm1
429            DO ji = fs_2, fs_jpim1   ! vector opt.
430               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
[1239]431            END DO
432         END DO
433      END DO
434
[6491]435      !                                 ! Save TKE prior to nn_etau addition 
436      e_niw(:,:,:) = en(:,:,:) 
437     
[1492]438      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
439      !                            !  TKE due to surface and internal wave breaking
440      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[6491]441      IF( nn_htau == 2 ) THEN           !* mixed-layer depth dependant length scale
442         DO jj = 2, jpjm1
443            DO ji = fs_2, fs_jpim1   ! vector opt.
444               htau(ji,jj) = rhtau * hmlp(ji,jj)
445            END DO
446         END DO
447      ENDIF
448#if defined key_iomput
449      !
450      CALL iom_put( "htau", htau(:,:) )  ! Check htau (even if constant in time)
451#endif
452      !
[2528]453      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
[1492]454         DO jk = 2, jpkm1
[1239]455            DO jj = 2, jpjm1
456               DO ji = fs_2, fs_jpim1   ! vector opt.
[1492]457                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
[5120]458                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1239]459               END DO
460            END DO
[1492]461         END DO
[2528]462      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
[1492]463         DO jj = 2, jpjm1
464            DO ji = fs_2, fs_jpim1   ! vector opt.
465               jk = nmln(ji,jj)
466               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
[5120]467                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1239]468            END DO
[1492]469         END DO
[2528]470      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
[1705]471!CDIR NOVERRCHK
472         DO jk = 2, jpkm1
473!CDIR NOVERRCHK
474            DO jj = 2, jpjm1
475!CDIR NOVERRCHK
476               DO ji = fs_2, fs_jpim1   ! vector opt.
477                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
478                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
[4990]479                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
[2528]480                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
481                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
[1705]482                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
[5120]483                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1705]484               END DO
485            END DO
486         END DO
[6491]487      ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up)
488         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau
489            DO jj = 2, jpjm1
490               DO ji = fs_2, fs_jpim1   ! vector opt.
491                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) )
492               END DO
493            END DO
494         ENDIF
495         DO jk = 2, jpkm1
496            DO jj = 2, jpjm1
497               DO ji = fs_2, fs_jpim1   ! vector opt.
498                  en(ji,jj,jk) = en(ji,jj,jk) + efr(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
499                     &                                                   * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
500               END DO
501            END DO
502         END DO
[1239]503      ENDIF
[1492]504      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
505      !
[6491]506      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term 
507         DO jj = 2, jpjm1 
508            DO ji = fs_2, fs_jpim1   ! vector opt. 
509               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 
510            END DO 
511         END DO 
512      END DO 
513     
514      CALL lbc_lnk( e_niw, 'W', 1. ) 
515      !
[3294]516      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
517      CALL wrk_dealloc( jpi,jpj, zhlc ) 
518      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
[2715]519      !
[3294]520      IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
521      !
[1239]522   END SUBROUTINE tke_tke
523
[1492]524
525   SUBROUTINE tke_avn
[1239]526      !!----------------------------------------------------------------------
[1492]527      !!                   ***  ROUTINE tke_avn  ***
[1239]528      !!
[1492]529      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
530      !!
531      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
532      !!              the tke_tke routine). First, the now mixing lenth is
533      !!      computed from en and the strafification (N^2), then the mixings
534      !!      coefficients are computed.
535      !!              - Mixing length : a first evaluation of the mixing lengh
536      !!      scales is:
537      !!                      mxl = sqrt(2*en) / N 
538      !!      where N is the brunt-vaisala frequency, with a minimum value set
[2528]539      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
[1492]540      !!        The mixing and dissipative length scale are bound as follow :
541      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
542      !!                        zmxld = zmxlm = mxl
543      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
544      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
545      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
546      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
547      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
548      !!                    the surface to obtain ldown. the resulting length
549      !!                    scales are:
550      !!                        zmxld = sqrt( lup * ldown )
551      !!                        zmxlm = min ( lup , ldown )
552      !!              - Vertical eddy viscosity and diffusivity:
553      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
554      !!                      avt = max( avmb, pdlr * avm ) 
555      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
556      !!
557      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
558      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
[1239]559      !!----------------------------------------------------------------------
[2715]560      INTEGER  ::   ji, jj, jk   ! dummy loop indices
561      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
562      REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      -
563      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
[3294]564      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
[1239]565      !!--------------------------------------------------------------------
[3294]566      !
567      IF( nn_timing == 1 )  CALL timing_start('tke_avn')
[1239]568
[3294]569      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
570
[1492]571      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
572      !                     !  Mixing length
573      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
574      !
575      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
576      !
[5120]577      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
578      zmxlm(:,:,:)  = rmxl_min   
579      zmxld(:,:,:)  = rmxl_min
580      !
[2528]581      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
[4990]582         DO jj = 2, jpjm1
583            DO ji = fs_2, fs_jpim1
[5120]584               zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
585               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
[4990]586            END DO
587         END DO
588      ELSE
[5120]589         zmxlm(:,:,1) = rn_mxl0
[1239]590      ENDIF
591      !
592!CDIR NOVERRCHK
[5120]593      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
[1239]594!CDIR NOVERRCHK
[5120]595         DO jj = 2, jpjm1
596!CDIR NOVERRCHK
597            DO ji = fs_2, fs_jpim1   ! vector opt.
[1239]598               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
[5120]599               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
[1239]600            END DO
601         END DO
602      END DO
[1492]603      !
604      !                     !* Physical limits for the mixing length
605      !
[5120]606      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value
[2528]607      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
[1492]608      !
[1239]609      SELECT CASE ( nn_mxl )
610      !
[5120]611      ! where wmask = 0 set zmxlm == fse3w
[1239]612      CASE ( 0 )           ! bounded by the distance to surface and bottom
[5120]613         DO jk = 2, jpkm1
614            DO jj = 2, jpjm1
615               DO ji = fs_2, fs_jpim1   ! vector opt.
[4990]616                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
[2528]617                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
[5120]618                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
619                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
620                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
[1239]621               END DO
622            END DO
623         END DO
624         !
625      CASE ( 1 )           ! bounded by the vertical scale factor
[5120]626         DO jk = 2, jpkm1
627            DO jj = 2, jpjm1
628               DO ji = fs_2, fs_jpim1   ! vector opt.
[1239]629                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
630                  zmxlm(ji,jj,jk) = zemxl
631                  zmxld(ji,jj,jk) = zemxl
632               END DO
633            END DO
634         END DO
635         !
636      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
[5120]637         DO jk = 2, jpkm1         ! from the surface to the bottom :
638            DO jj = 2, jpjm1
639               DO ji = fs_2, fs_jpim1   ! vector opt.
[1239]640                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
641               END DO
[5120]642            END DO
643         END DO
644         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
645            DO jj = 2, jpjm1
646               DO ji = fs_2, fs_jpim1   ! vector opt.
[1239]647                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
648                  zmxlm(ji,jj,jk) = zemxl
649                  zmxld(ji,jj,jk) = zemxl
650               END DO
651            END DO
652         END DO
653         !
654      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
[5120]655         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
656            DO jj = 2, jpjm1
657               DO ji = fs_2, fs_jpim1   ! vector opt.
[1239]658                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
659               END DO
[5120]660            END DO
661         END DO
662         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
663            DO jj = 2, jpjm1
664               DO ji = fs_2, fs_jpim1   ! vector opt.
[1239]665                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
666               END DO
667            END DO
668         END DO
669!CDIR NOVERRCHK
670         DO jk = 2, jpkm1
671!CDIR NOVERRCHK
672            DO jj = 2, jpjm1
673!CDIR NOVERRCHK
674               DO ji = fs_2, fs_jpim1   ! vector opt.
675                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
676                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
677                  zmxlm(ji,jj,jk) = zemlm
678                  zmxld(ji,jj,jk) = zemlp
679               END DO
680            END DO
681         END DO
682         !
683      END SELECT
[1492]684      !
[1239]685# if defined key_c1d
[1492]686      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
[1239]687      e_mix(:,:,:) = zmxlm(:,:,:)
688# endif
689
[1492]690      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
691      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
692      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1239]693!CDIR NOVERRCHK
[1492]694      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
[1239]695!CDIR NOVERRCHK
696         DO jj = 2, jpjm1
697!CDIR NOVERRCHK
698            DO ji = fs_2, fs_jpim1   ! vector opt.
699               zsqen = SQRT( en(ji,jj,jk) )
700               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
[5120]701               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
702               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
[1239]703               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
704            END DO
705         END DO
706      END DO
[1492]707      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
708      !
[5120]709      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
710         DO jj = 2, jpjm1
711            DO ji = fs_2, fs_jpim1   ! vector opt.
712               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
713               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
[4990]714            END DO
[1239]715         END DO
716      END DO
717      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
[1492]718      !
719      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
[5120]720         DO jk = 2, jpkm1
721            DO jj = 2, jpjm1
722               DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]723                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
[1492]724                  !                                          ! shear
725                  zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   &
726                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
727                  zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   &
728                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
729                  !                                          ! local Richardson number
[2528]730                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
731                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
[1492]732!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
[2528]733!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
[5120]734                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
[1492]735# if defined key_c1d
[5120]736                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number
737                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri
[1239]738# endif
739              END DO
740            END DO
741         END DO
742      ENDIF
743      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
744
745      IF(ln_ctl) THEN
746         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
747         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
748            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
749      ENDIF
750      !
[3294]751      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
752      !
753      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
754      !
[1492]755   END SUBROUTINE tke_avn
[1239]756
[1492]757
[2528]758   SUBROUTINE zdf_tke_init
[1239]759      !!----------------------------------------------------------------------
[2528]760      !!                  ***  ROUTINE zdf_tke_init  ***
[1239]761      !!                     
762      !! ** Purpose :   Initialization of the vertical eddy diffivity and
[1492]763      !!              viscosity when using a tke turbulent closure scheme
[1239]764      !!
[1601]765      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
[1492]766      !!              called at the first timestep (nit000)
[1239]767      !!
[1601]768      !! ** input   :   Namlist namzdf_tke
[1239]769      !!
770      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
771      !!----------------------------------------------------------------------
772      INTEGER ::   ji, jj, jk   ! dummy loop indices
[6487]773      INTEGER ::   ios, ierr
[1239]774      !!
[2528]775      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
776         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
777         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
[6491]778         &                 nn_etau , nn_htau  , rn_efr , rn_c   
[1239]779      !!----------------------------------------------------------------------
[6491]780
[4147]781      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
782      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
783901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
784
785      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
786      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
787902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
[10759]788      IF(lwm .AND. nprint > 2) WRITE ( numond, namzdf_tke )
[2715]789      !
[2528]790      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
[6498]791# if defined key_zdftmx_new
792      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
793      rn_emin  = 1.e-10_wp
794      rmxl_min = 1.e-03_wp
795      IF(lwp) THEN                  ! Control print
796         WRITE(numout,*)
797         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
798         WRITE(numout,*) '~~~~~~~~~~~~'
[10774]799         IF(lflush) CALL flush(numout)
[6498]800      ENDIF
801# else
[2528]802      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
[6498]803# endif
[2715]804      !
[1492]805      IF(lwp) THEN                    !* Control print
[1239]806         WRITE(numout,*)
[2528]807         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
808         WRITE(numout,*) '~~~~~~~~~~~~'
[1601]809         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
[1705]810         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
811         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
812         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
813         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
814         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
815         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
816         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
817         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
818         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
[2528]819         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
820         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
821         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
[1705]822         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
823         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
824         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
[6491]825         WRITE(numout,*) '      fraction of TKE added within the mixed layer by nn_etau rn_c    = ', rn_c
[1239]826         WRITE(numout,*)
[1601]827         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
[10774]828         IF(lflush) CALL flush(numout)
[1239]829      ENDIF
[2715]830      !
831      !                              ! allocate tke arrays
832      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
833      !
[1492]834      !                               !* Check of some namelist values
[4990]835      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
836      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
[6491]837      IF( nn_htau < 0  .OR.  nn_htau > 5 )   CALL ctl_stop( 'bad flag: nn_htau is 0 to 5    ' )
[5407]838      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
[1239]839
[2528]840      IF( ln_mxl0 ) THEN
841         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
[10774]842         IF(lwp .AND. lflush) CALL flush(numout)
[2528]843         rn_mxl0 = rmxl_min
844      ENDIF
845     
[6487]846      IF( nn_etau == 2  ) THEN
847          ierr = zdf_mxl_alloc()
848          nmln(:,:) = nlb10           ! Initialization of nmln
849      ENDIF
[1239]850
[6491]851      IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN
852          ierr = zdf_mxl_alloc()
853          nmln(:,:) = nlb10           ! Initialization of nmln
854      ENDIF
855
[1492]856      !                               !* depth of penetration of surface tke
857      IF( nn_etau /= 0 ) THEN     
[6491]858         htau(:,:) = 0._wp
[1601]859         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
[2528]860         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
861            htau(:,:) = 10._wp
862         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
863            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
[6491]864         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer
865            rhtau = -1._wp / LOG( 1._wp - rn_c )
866         CASE( 3 )                                 ! F(latitude) : 0.5m to 15m poleward of 20 degrees
867            htau(:,:) = MAX(  0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
868         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south
869            DO jj = 2, jpjm1
870               DO ji = fs_2, fs_jpim1   ! vector opt.
871                  IF( gphit(ji,jj) <= 0._wp ) THEN
872                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
873                  ELSE
874                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
875                  ENDIF
876               END DO
877            END DO
878         CASE ( 5 )                                ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south,
879            DO jj = 2, jpjm1                       !               10m to 30m between 30/45 degrees south
880               DO ji = fs_2, fs_jpim1   ! vector opt.
881                  IF( gphit(ji,jj) <= -30._wp ) THEN
882                     htau(ji,jj) = MAX(  10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) )   )
883                  ELSE
884                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
885                  ENDIF
886               END DO
887            END DO
[1492]888         END SELECT
[6491]889         !
890         IF( nn_etau == 4 .AND. nn_htau /= 2 ) THEN            ! efr dependant on constant htau
891            DO jj = 2, jpjm1
892               DO ji = fs_2, fs_jpim1   ! vector opt.
893                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) )
894               END DO
895            END DO
896         ENDIF
[1492]897      ENDIF
898      !                               !* set vertical eddy coef. to the background value
[1239]899      DO jk = 1, jpk
[5120]900         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
901         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
902         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
903         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
[1239]904      END DO
[2528]905      dissl(:,:,:) = 1.e-12_wp
[2715]906      !                             
907      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
[1239]908      !
[2528]909   END SUBROUTINE zdf_tke_init
[1239]910
911
[1531]912   SUBROUTINE tke_rst( kt, cdrw )
[1239]913     !!---------------------------------------------------------------------
[1531]914     !!                   ***  ROUTINE tke_rst  ***
[1239]915     !!                     
916     !! ** Purpose :   Read or write TKE file (en) in restart file
917     !!
918     !! ** Method  :   use of IOM library
919     !!                if the restart does not contain TKE, en is either
[1537]920     !!                set to rn_emin or recomputed
[1239]921     !!----------------------------------------------------------------------
[2715]922     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
923     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
[1239]924     !
[1481]925     INTEGER ::   jit, jk   ! dummy loop indices
[2715]926     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
[1239]927     !!----------------------------------------------------------------------
928     !
[1481]929     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
930        !                                   ! ---------------
931        IF( ln_rstart ) THEN                   !* Read the restart file
932           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
933           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
934           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
935           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
936           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
937           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
938           !
939           IF( id1 > 0 ) THEN                       ! 'en' exists
[1239]940              CALL iom_get( numror, jpdom_autoglo, 'en', en )
[1481]941              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
[9321]942                 IF(nn_timing == 2)  CALL timing_start('iom_rstget')
[1481]943                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
944                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
945                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
946                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
947                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
[9321]948                 IF(nn_timing == 2)  CALL timing_stop('iom_rstget')
[1492]949              ELSE                                                 ! one at least array is missing
950                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
[1481]951              ENDIF
952           ELSE                                     ! No TKE array found: initialisation
953              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
[10774]954              IF(lwp .AND. lflush) CALL flush(numout)
[1239]955              en (:,:,:) = rn_emin * tmask(:,:,:)
[1492]956              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
[5112]957              !
958              avt_k (:,:,:) = avt (:,:,:)
959              avm_k (:,:,:) = avm (:,:,:)
960              avmu_k(:,:,:) = avmu(:,:,:)
961              avmv_k(:,:,:) = avmv(:,:,:)
962              !
[1531]963              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
[1239]964           ENDIF
[1481]965        ELSE                                   !* Start from rest
966           en(:,:,:) = rn_emin * tmask(:,:,:)
967           DO jk = 1, jpk                           ! set the Kz to the background value
[5120]968              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
969              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
970              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
971              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
[1481]972           END DO
[1239]973        ENDIF
[1481]974        !
975     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
976        !                                   ! -------------------
[10774]977        IF(lwp .AND. nprint > 0) THEN
978           WRITE(numout,*) '---- tke-rst ----'
979           IF(lflush) CALL flush(numout)
980        ENDIF
[9321]981        IF(nn_timing == 2)  CALL timing_start('iom_rstput')
[3632]982        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
983        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
984        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
985        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
986        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
987        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
[9321]988        IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
[1481]989        !
[1239]990     ENDIF
991     !
[1531]992   END SUBROUTINE tke_rst
[1239]993
994#else
995   !!----------------------------------------------------------------------
996   !!   Dummy module :                                        NO TKE scheme
997   !!----------------------------------------------------------------------
[1531]998   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
[1239]999CONTAINS
[2528]1000   SUBROUTINE zdf_tke_init           ! Dummy routine
1001   END SUBROUTINE zdf_tke_init
1002   SUBROUTINE zdf_tke( kt )          ! Dummy routine
[1531]1003      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
1004   END SUBROUTINE zdf_tke
1005   SUBROUTINE tke_rst( kt, cdrw )
[1492]1006     CHARACTER(len=*) ::   cdrw
[1531]1007     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
1008   END SUBROUTINE tke_rst
[1239]1009#endif
1010
1011   !!======================================================================
[1531]1012END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.