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

source: branches/2012/dev_LOCEAN_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 3653

Last change on this file since 3653 was 3653, checked in by cetlod, 11 years ago

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