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

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 5098

Last change on this file since 5098 was 5098, checked in by mathiot, 9 years ago

add wmask, wumask, wvmask and restore loop order and add flag to ignore specific isf code if no isf

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