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

source: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 4644

Last change on this file since 4644 was 4644, checked in by acc, 10 years ago

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