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/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 4861

Last change on this file since 4861 was 4861, checked in by rblod, 9 years ago

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