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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 9383

Last change on this file since 9383 was 9383, checked in by andmirek, 6 years ago

#2050 fixes and changes

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