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

source: branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 8055

Last change on this file since 8055 was 8055, checked in by gm, 7 years ago

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