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

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 8519

Last change on this file since 8519 was 8519, checked in by lovato, 7 years ago

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