New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdftke.F90 in branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 3524

Last change on this file since 3524 was 3524, checked in by gm, 11 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. add USE lib_fortran when SIGN is used (TOP,OPA,LIM2&3) ; salt flux names start with sfx_ in LIM3

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