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

source: branches/UKMO/test_moci_test_suite/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 8243

Last change on this file since 8243 was 8243, checked in by andmirek, 7 years ago

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