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

Last change on this file since 3186 was 3186, checked in by smasson, 10 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo

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