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

source: branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 11384

Last change on this file since 11384 was 11384, checked in by mattmartin, 5 years ago

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