source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 3229

Last change on this file since 3229 was 3229, checked in by charris, 10 years ago

Added timing calls to most significant routines in LDF, SBC and ZDF.

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