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

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 8863

Last change on this file since 8863 was 8863, checked in by flavoni, 6 years ago

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