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

source: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 1617

Last change on this file since 1617 was 1617, checked in by ctlod, 15 years ago

remove the TKE profile of penetration case (2) for the parameter nn_htau, see ticket: #521
File size: 40.4 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   !!----------------------------------------------------------------------
28#if defined key_zdftke   ||   defined key_esopa
29   !!----------------------------------------------------------------------
30   !!   'key_zdftke'                                   TKE vertical physics
31   !!----------------------------------------------------------------------
32   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
33   !!   tke_tke       : tke time stepping: update tke at now time step (en)
34   !!   tke_avn       : compute mixing length scale and deduce avm and avt
35   !!   tke_init      : initialization, namelist read, and parameters control
36   !!   tke_rst       : read/write tke restart in ocean restart file
37   !!----------------------------------------------------------------------
38   USE oce            ! ocean dynamics and active tracers
39   USE dom_oce        ! ocean space and time domain
40   USE domvvl         ! ocean space and time domain : variable volume layer
41   USE zdf_oce        ! ocean vertical physics
42   USE sbc_oce        ! surface boundary condition: ocean
43   USE phycst         ! physical constants
44   USE zdfmxl         ! mixed layer
45   USE restart        ! only for lrst_oce
46   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
47   USE prtctl         ! Print control
48   USE in_out_manager ! I/O manager
49   USE iom            ! I/O manager library
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC   zdf_tke    ! routine called in step module
55   PUBLIC   tke_rst    ! routine called in step module
56
57   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
58
59#if defined key_c1d
60   !                                                                !!* 1D cfg only *
61   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales
62   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers
63#endif
64
65   !                                       !!! ** Namelist  namzdf_tke  **
66   LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not
67   INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3)
68   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length   [m]
69   REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length   [m]
70   INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1)
71   REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
72   REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation
73   REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke
74   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke           [m2/s2]
75   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke   [m2/s2]
76   REAL(wp) ::   rn_bshear= 1.e-20          ! background shear (>0)
77   INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2)
78   INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1)
79   REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean
80   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
81   REAL(wp) ::   rn_lc    = 0.15_wp         ! 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
85   REAL(wp), DIMENSION(jpi,jpj)     ::   htau      ! depth of tke penetration (nn_htau)
86   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   en        ! now turbulent kinetic energy   [m2/s2]
87   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dissl     ! now mixing lenght of dissipation
88
89   !! * Substitutions
90#  include "domzgr_substitute.h90"
91#  include "vectopt_loop_substitute.h90"
92   !!----------------------------------------------------------------------
93   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
94   !! $Id: zdftke2.F90 1201 2008-09-24 13:24:21Z rblod $
95   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
96   !!----------------------------------------------------------------------
97
98CONTAINS
99
100   SUBROUTINE zdf_tke( kt )
101      !!----------------------------------------------------------------------
102      !!                   ***  ROUTINE zdf_tke  ***
103      !!
104      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
105      !!              coefficients using a turbulent closure scheme (TKE).
106      !!
107      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
108      !!              is computed from a prognostic equation :
109      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
110      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
111      !!                  + avt N^2                      ! stratif. destruc.
112      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
113      !!      with the boundary conditions:
114      !!         surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) )
115      !!         bottom : en = rn_emin
116      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
117      !!
118      !!        The now Turbulent kinetic energy is computed using the following
119      !!      time stepping: implicit for vertical diffusion term, linearized semi
120      !!      implicit for kolmogoroff dissipation term, and explicit forward for
121      !!      both buoyancy and shear production terms. Therefore a tridiagonal
122      !!      linear system is solved. Note that buoyancy and shear terms are
123      !!      discretized in a energy conserving form (Bruchard 2002).
124      !!
125      !!        The dissipative and mixing length scale are computed from en and
126      !!      the stratification (see tke_avn)
127      !!
128      !!        The now vertical eddy vicosity and diffusivity coefficients are
129      !!      given by:
130      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
131      !!              avt = max( avmb, pdl * avm                 ) 
132      !!              eav = max( avmb, avm )
133      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
134      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
135      !!
136      !! ** Action  :   compute en (now turbulent kinetic energy)
137      !!                update avt, avmu, avmv (before vertical eddy coef.)
138      !!
139      !! References : Gaspar et al., JGR, 1990,
140      !!              Blanke and Delecluse, JPO, 1991
141      !!              Mellor and Blumberg, JPO 2004
142      !!              Axell, JGR, 2002
143      !!              Bruchard OM 2002
144      !!----------------------------------------------------------------------
145      INTEGER, INTENT(in) ::   kt   ! ocean time step
146      !!----------------------------------------------------------------------
147      !
148      IF( kt == nit000 )   CALL tke_init     ! initialisation
149                           !
150                           CALL tke_tke      ! now tke (en)
151                           !
152                           CALL tke_avn      ! now avt, avm, avmu, avmv
153      !
154   END SUBROUTINE zdf_tke
155
156
157   SUBROUTINE tke_tke
158      !!----------------------------------------------------------------------
159      !!                   ***  ROUTINE tke_tke  ***
160      !!
161      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
162      !!
163      !! ** Method  : - TKE surface boundary condition
164      !!              - source term due to Langmuir cells (ln_lc=T)
165      !!              - source term due to shear (saved in avmu, avmv arrays)
166      !!              - Now TKE : resolution of the TKE equation by inverting
167      !!                a tridiagonal linear system by a "methode de chasse"
168      !!              - increase TKE due to surface and internal wave breaking
169      !!
170      !! ** Action  : - en : now turbulent kinetic energy)
171      !!              - avmu, avmv : production of TKE by shear at u and v-points
172      !!                (= Kz dz[Ub] * dz[Un] )
173      !! ---------------------------------------------------------------------
174      USE oce,   zdiag  =>   ua   ! use ua as workspace
175      USE oce,   zd_up  =>   va   ! use va as workspace
176      USE oce,   zd_lw  =>   ta   ! use ta as workspace
177      !!
178      INTEGER  ::   ji, jj, jk                ! dummy loop arguments
179      REAL(wp) ::   zbbrau, zesurf, zesh2     ! temporary scalars
180      REAL(wp) ::   zfact1, zfact2, zfact3    !    -         -
181      REAL(wp) ::   ztx2  , zty2  , zcof      !    -         -
182      REAL(wp) ::   zus   , zwlc  , zind      !    -         -
183      REAL(wp) ::   zzd_up, zzd_lw            !    -         -
184      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
185      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
186      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
187      !!--------------------------------------------------------------------
188      !
189      zbbrau =  .5 * rn_ebb / rau0       ! Local constant initialisation
190      zfact1 = -.5 * rdt 
191      zfact2 = 1.5 * rdt * rn_ediss
192      zfact3 = 0.5       * rn_ediss
193      !
194      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
195      !                     !  Surface boundary condition on tke
196      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
197!CDIR NOVERRCHK
198      DO jj = 2, jpjm1         ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
199!CDIR NOVERRCHK
200         DO ji = fs_2, fs_jpim1   ! vector opt.
201            ztx2   = utau(ji-1,jj  ) + utau(ji,jj)
202            zty2   = vtau(ji  ,jj-1) + vtau(ji,jj)
203            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
204            en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
205         END DO
206      END DO
207      !
208      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
209      !                     !  Bottom boundary condition on tke
210      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
211!!gm to be added soon
212      !
213      !
214      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
215      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke
216         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
217         !
218         !                        !* total energy produce by LC : cumulative sum over jk
219         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
220         DO jk = 2, jpk
221            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
222         END DO
223         !                        !* finite Langmuir Circulation depth
224         imlc(:,:) = mbathy(:,:)         ! Initialization to the number of w ocean point mbathy
225         DO jk = jpkm1, 2, -1
226            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
227               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
228                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
229                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
230               END DO
231            END DO
232         END DO
233         !                               ! finite LC depth
234# if defined key_vectopt_loop
235         DO jj = 1, 1
236            DO ji = 1, jpij   ! vector opt. (forced unrolling)
237# else
238         DO jj = 1, jpj 
239            DO ji = 1, jpi
240# endif
241               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
242            END DO
243         END DO
244# if defined key_c1d
245         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
246# endif
247!CDIR NOVERRCHK
248         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
249!CDIR NOVERRCHK
250            DO jj = 2, jpjm1
251!CDIR NOVERRCHK
252               DO ji = fs_2, fs_jpim1   ! vector opt.
253!!gm replace here wndn by a formulation with the stress module
254                  zus  = 0.016 * wndm(ji,jj)                  ! Stokes drift
255                  !                                           ! vertical velocity due to LC
256                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
257                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
258                  !                                           ! TKE Langmuir circulation source term
259                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
260               END DO
261            END DO
262         END DO
263         !
264      ENDIF
265      !
266      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
267      !                     !  Now Turbulent kinetic energy (output in en)
268      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
269      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
270      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
271      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
272      !
273      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
274         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
275            DO ji = 1, jpi
276               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
277                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
278                  &           / (  fse3uw_n(ji,jj,jk)         &
279                  &              * fse3uw_b(ji,jj,jk) )
280               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
281                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
282                  &                            / (  fse3vw_n(ji,jj,jk)               &
283                  &                              *  fse3vw_b(ji,jj,jk)  )
284            END DO
285         END DO
286      END DO
287      !
288      DO jk = 2, jpkm1           !* Matrix and right hand side in en
289         DO jj = 2, jpjm1
290            DO ji = fs_2, fs_jpim1   ! vector opt.
291               zcof   = zfact1 * tmask(ji,jj,jk)
292               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
293                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
294               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
295                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
296                  !                                                           ! shear prod. at w-point weightened by mask
297               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
298                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
299                  !
300               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
301               zd_lw(ji,jj,jk) = zzd_lw
302               zdiag(ji,jj,jk) = 1.e0 - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
303               !
304               !                                   ! right hand side in en
305               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
306                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk)
307            END DO
308         END DO
309      END DO
310      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
311      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
312         DO jj = 2, jpjm1
313            DO ji = fs_2, fs_jpim1    ! vector opt.
314               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
315            END DO
316         END DO
317      END DO
318      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
319         DO ji = fs_2, fs_jpim1    ! vector opt.
320            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
321         END DO
322      END DO
323      DO jk = 3, jpkm1
324         DO jj = 2, jpjm1
325            DO ji = fs_2, fs_jpim1    ! vector opt.
326               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)
327            END DO
328         END DO
329      END DO
330      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
331         DO ji = fs_2, fs_jpim1    ! vector opt.
332            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
333         END DO
334      END DO
335      DO jk = jpk-2, 2, -1
336         DO jj = 2, jpjm1
337            DO ji = fs_2, fs_jpim1    ! vector opt.
338               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
339            END DO
340         END DO
341      END DO
342      DO jk = 2, jpkm1                             ! set the minimum value of tke
343         DO jj = 2, jpjm1
344            DO ji = fs_2, fs_jpim1   ! vector opt.
345               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
346            END DO
347         END DO
348      END DO
349
350      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
351      !                            !  TKE due to surface and internal wave breaking
352      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
353      IF( nn_etau == 1 ) THEN           !*  penetration throughout the water column
354         DO jk = 2, jpkm1
355            DO jj = 2, jpjm1
356               DO ji = fs_2, fs_jpim1   ! vector opt.
357                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
358                     &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
359               END DO
360            END DO
361         END DO
362      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)
363         DO jj = 2, jpjm1
364            DO ji = fs_2, fs_jpim1   ! vector opt.
365               jk = nmln(ji,jj)
366               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
367                  &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
368            END DO
369         END DO
370      ENDIF
371      !
372      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
373      !
374   END SUBROUTINE tke_tke
375
376
377   SUBROUTINE tke_avn
378      !!----------------------------------------------------------------------
379      !!                   ***  ROUTINE tke_avn  ***
380      !!
381      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
382      !!
383      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
384      !!              the tke_tke routine). First, the now mixing lenth is
385      !!      computed from en and the strafification (N^2), then the mixings
386      !!      coefficients are computed.
387      !!              - Mixing length : a first evaluation of the mixing lengh
388      !!      scales is:
389      !!                      mxl = sqrt(2*en) / N 
390      !!      where N is the brunt-vaisala frequency, with a minimum value set
391      !!      to rn_lmin (rn_lmin0) in the interior (surface) ocean.
392      !!        The mixing and dissipative length scale are bound as follow :
393      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
394      !!                        zmxld = zmxlm = mxl
395      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
396      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
397      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
398      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
399      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
400      !!                    the surface to obtain ldown. the resulting length
401      !!                    scales are:
402      !!                        zmxld = sqrt( lup * ldown )
403      !!                        zmxlm = min ( lup , ldown )
404      !!              - Vertical eddy viscosity and diffusivity:
405      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
406      !!                      avt = max( avmb, pdlr * avm ) 
407      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
408      !!
409      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
410      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
411      !!----------------------------------------------------------------------
412      USE oce,     zmpdl  =>   ua   ! use ua as workspace
413      USE oce,     zmxlm  =>   va   ! use va as workspace
414      USE oce,     zmxld  =>   ta   ! use ta as workspace
415      !!
416      INTEGER  ::   ji, jj, jk            ! dummy loop arguments
417      REAL(wp) ::   zrn2, zraug           ! temporary scalars
418      REAL(wp) ::   ztx2, zdku            !    -         -
419      REAL(wp) ::   zty2, zdkv            !    -         -
420      REAL(wp) ::   zcoef, zav            !    -         -
421      REAL(wp) ::   zpdlr, zri, zsqen     !    -         -
422      REAL(wp) ::   zemxl, zemlm, zemlp   !    -         -
423      !!--------------------------------------------------------------------
424
425      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
426      !                     !  Mixing length
427      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
428      !
429      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
430      !
431      IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=vkarmn*2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
432!!gm  this should be useless
433         zmxlm(:,:,1) = 0.e0
434!!gm end
435         zraug = 0.5 * vkarmn * 2.e5 / ( rau0 * grav )
436         DO jj = 2, jpjm1
437!CDIR NOVERRCHK
438            DO ji = fs_2, fs_jpim1   ! vector opt.
439               ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
440               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
441               zmxlm(ji,jj,1) = MAX(  rn_lmin0,  zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )  )
442            END DO
443         END DO
444      ELSE                       ! surface set to the minimum value
445         zmxlm(:,:,1) = rn_lmin0
446      ENDIF
447      zmxlm(:,:,jpk) = rn_lmin   ! bottom set to the interior minium value
448      !
449!CDIR NOVERRCHK
450      DO jk = 2, jpkm1           ! interior value : l=sqrt(2*e/n**2)
451!CDIR NOVERRCHK
452         DO jj = 2, jpjm1
453!CDIR NOVERRCHK
454            DO ji = fs_2, fs_jpim1   ! vector opt.
455               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
456               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
457            END DO
458         END DO
459      END DO
460      !
461!!gm  CAUTION: to be added here the bottom boundary condition on zmxlm
462      !
463      !                     !* Physical limits for the mixing length
464      !
465      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
466      zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value
467      !
468      SELECT CASE ( nn_mxl )
469      !
470      CASE ( 0 )           ! bounded by the distance to surface and bottom
471         DO jk = 2, jpkm1
472            DO jj = 2, jpjm1
473               DO ji = fs_2, fs_jpim1   ! vector opt.
474                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
475                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
476                  zmxlm(ji,jj,jk) = zemxl
477                  zmxld(ji,jj,jk) = zemxl
478               END DO
479            END DO
480         END DO
481         !
482      CASE ( 1 )           ! bounded by the vertical scale factor
483         DO jk = 2, jpkm1
484            DO jj = 2, jpjm1
485               DO ji = fs_2, fs_jpim1   ! vector opt.
486                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
487                  zmxlm(ji,jj,jk) = zemxl
488                  zmxld(ji,jj,jk) = zemxl
489               END DO
490            END DO
491         END DO
492         !
493      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
494         DO jk = 2, jpkm1         ! from the surface to the bottom :
495            DO jj = 2, jpjm1
496               DO ji = fs_2, fs_jpim1   ! vector opt.
497                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
498               END DO
499            END DO
500         END DO
501         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
502            DO jj = 2, jpjm1
503               DO ji = fs_2, fs_jpim1   ! vector opt.
504                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
505                  zmxlm(ji,jj,jk) = zemxl
506                  zmxld(ji,jj,jk) = zemxl
507               END DO
508            END DO
509         END DO
510         !
511      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
512         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
513            DO jj = 2, jpjm1
514               DO ji = fs_2, fs_jpim1   ! vector opt.
515                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
516               END DO
517            END DO
518         END DO
519         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
520            DO jj = 2, jpjm1
521               DO ji = fs_2, fs_jpim1   ! vector opt.
522                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
523               END DO
524            END DO
525         END DO
526!CDIR NOVERRCHK
527         DO jk = 2, jpkm1
528!CDIR NOVERRCHK
529            DO jj = 2, jpjm1
530!CDIR NOVERRCHK
531               DO ji = fs_2, fs_jpim1   ! vector opt.
532                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
533                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
534                  zmxlm(ji,jj,jk) = zemlm
535                  zmxld(ji,jj,jk) = zemlp
536               END DO
537            END DO
538         END DO
539         !
540      END SELECT
541      !
542# if defined key_c1d
543      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
544      e_mix(:,:,:) = zmxlm(:,:,:)
545# endif
546
547      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
548      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
549      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
550!CDIR NOVERRCHK
551      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
552!CDIR NOVERRCHK
553         DO jj = 2, jpjm1
554!CDIR NOVERRCHK
555            DO ji = fs_2, fs_jpim1   ! vector opt.
556               zsqen = SQRT( en(ji,jj,jk) )
557               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
558               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
559               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
560               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
561            END DO
562         END DO
563      END DO
564      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
565      !
566      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
567         DO jj = 2, jpjm1
568            DO ji = fs_2, fs_jpim1   ! vector opt.
569               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
570               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
571            END DO
572         END DO
573      END DO
574      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
575      !
576      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
577         DO jk = 2, jpkm1
578            DO jj = 2, jpjm1
579               DO ji = fs_2, fs_jpim1   ! vector opt.
580                  zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
581                  !                                          ! shear
582                  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) )   &
583                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
584                  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) )   &
585                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
586                  !                                          ! local Richardson number
587                  zri   = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk) / (zcoef * (zdku + zdkv + rn_bshear ) )
588                  zpdlr = MAX(  0.1,  0.2 / MAX( 0.2 , zri )  )
589!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
590!!gm              zpdlr = MAX(  0.1,  ri_crit / MAX( ri_crit , zri )  )
591                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
592# if defined key_c1d
593                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
594                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
595# endif
596              END DO
597            END DO
598         END DO
599      ENDIF
600      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
601
602      IF(ln_ctl) THEN
603         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
604         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
605            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
606      ENDIF
607      !
608   END SUBROUTINE tke_avn
609
610
611   SUBROUTINE tke_init
612      !!----------------------------------------------------------------------
613      !!                  ***  ROUTINE tke_init  ***
614      !!                     
615      !! ** Purpose :   Initialization of the vertical eddy diffivity and
616      !!              viscosity when using a tke turbulent closure scheme
617      !!
618      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
619      !!              called at the first timestep (nit000)
620      !!
621      !! ** input   :   Namlist namzdf_tke
622      !!
623      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
624      !!----------------------------------------------------------------------
625      INTEGER ::   ji, jj, jk   ! dummy loop indices
626      !!
627      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb, rn_emin,   &
628         &                 rn_emin0, rn_bshear, nn_mxl, ln_mxl0,   &
629         &                 rn_lmin , rn_lmin0 , nn_pdl, nn_etau,   &
630         &                 nn_htau , rn_efr   , ln_lc , rn_lc 
631      !!----------------------------------------------------------------------
632
633      REWIND ( numnam )               !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
634      READ   ( numnam, namzdf_tke )
635     
636      ri_cri = 2. / ( 2. + rn_ediss / rn_ediff )      ! resulting critical Richardson number
637
638      IF(lwp) THEN                    !* Control print
639         WRITE(numout,*)
640         WRITE(numout,*) 'zdf_tke : tke turbulent closure scheme - initialisation'
641         WRITE(numout,*) '~~~~~~~~'
642         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
643         WRITE(numout,*) '      coef. to compute avt                        rn_ediff = ', rn_ediff
644         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss = ', rn_ediss
645         WRITE(numout,*) '      tke surface input coef.                     rn_ebb   = ', rn_ebb
646         WRITE(numout,*) '      minimum value of tke                        rn_emin  = ', rn_emin
647         WRITE(numout,*) '      surface minimum value of tke                rn_emin0 = ', rn_emin0
648         WRITE(numout,*) '      background shear (>0)                       rn_bshear= ', rn_bshear
649         WRITE(numout,*) '      mixing length type                          nn_mxl   = ', nn_mxl
650         WRITE(numout,*) '      prandl number flag                          nn_pdl   = ', nn_pdl
651         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0  = ', ln_mxl0
652         WRITE(numout,*) '      surface  mixing length minimum value        rn_lmin0 = ', rn_lmin0
653         WRITE(numout,*) '      interior mixing length minimum value        rn_lmin0 = ', rn_lmin0
654         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau  = ', nn_etau
655         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau  = ', nn_htau
656         WRITE(numout,*) '      % of rn_emin0 which pene. the thermocline   rn_efr   = ', rn_efr
657         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc    = ', ln_lc
658         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc    = ', rn_lc
659         WRITE(numout,*)
660         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
661      ENDIF
662
663      !                               !* Check of some namelist values
664      IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
665      IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
666      IF( nn_htau < 0    .OR. nn_htau > 1   )   CALL ctl_stop( 'bad flag: nn_htau is 0 or 1    ' )
667      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' )
668
669      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
670
671      !                               !* depth of penetration of surface tke
672      IF( nn_etau /= 0 ) THEN     
673         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
674         CASE( 0 )                                    ! constant depth penetration (here 10 meters)
675            htau(:,:) = 10.e0
676         CASE( 1 )                                    ! F(latitude) : 0.5m to 30m at high lat.
677            DO jj = 1, jpj
678               DO ji = 1, jpi 
679                  htau(ji,jj) = MAX( 0.5, 3./4. * MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
680               END DO
681            END DO
682         END SELECT
683      ENDIF
684
685      !                               !* set vertical eddy coef. to the background value
686      DO jk = 1, jpk
687         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
688         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
689         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
690         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
691      END DO
692      dissl(:,:,:) = 1.e-12
693      !                               !* read or initialize all required files
694      CALL tke_rst( nit000, 'READ' )
695      !
696   END SUBROUTINE tke_init
697
698
699   SUBROUTINE tke_rst( kt, cdrw )
700     !!---------------------------------------------------------------------
701     !!                   ***  ROUTINE tke_rst  ***
702     !!                     
703     !! ** Purpose :   Read or write TKE file (en) in restart file
704     !!
705     !! ** Method  :   use of IOM library
706     !!                if the restart does not contain TKE, en is either
707     !!                set to rn_emin or recomputed
708     !!----------------------------------------------------------------------
709     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
710     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
711     !
712     INTEGER ::   jit, jk   ! dummy loop indices
713     INTEGER ::   id1, id2, id3, id4, id5, id6
714     !!----------------------------------------------------------------------
715     !
716     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
717        !                                   ! ---------------
718        IF( ln_rstart ) THEN                   !* Read the restart file
719           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
720           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
721           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
722           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
723           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
724           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
725           !
726           IF( id1 > 0 ) THEN                       ! 'en' exists
727              CALL iom_get( numror, jpdom_autoglo, 'en', en )
728              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
729                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
730                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
731                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
732                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
733                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
734              ELSE                                                 ! one at least array is missing
735                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
736              ENDIF
737           ELSE                                     ! No TKE array found: initialisation
738              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
739              en (:,:,:) = rn_emin * tmask(:,:,:)
740              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
741              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
742           ENDIF
743        ELSE                                   !* Start from rest
744           en(:,:,:) = rn_emin * tmask(:,:,:)
745           DO jk = 1, jpk                           ! set the Kz to the background value
746              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
747              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
748              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
749              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
750           END DO
751        ENDIF
752        !
753     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
754        !                                   ! -------------------
755        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
756        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
757        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt   )
758        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm   )
759        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu  )
760        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  )
761        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
762        !
763     ENDIF
764     !
765   END SUBROUTINE tke_rst
766
767#else
768   !!----------------------------------------------------------------------
769   !!   Dummy module :                                        NO TKE scheme
770   !!----------------------------------------------------------------------
771   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
772CONTAINS
773   SUBROUTINE zdf_tke( kt )          ! Empty routine
774      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
775   END SUBROUTINE zdf_tke
776   SUBROUTINE tke_rst( kt, cdrw )
777     CHARACTER(len=*) ::   cdrw
778     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
779   END SUBROUTINE tke_rst
780#endif
781
782   !!======================================================================
783END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.