New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdftke.F90 in branches/UKMO/dev_merge_2017_restart_datestamp_GO6_mixing/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/dev_merge_2017_restart_datestamp_GO6_mixing/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 9497

Last change on this file since 9497 was 9497, checked in by davestorkey, 6 years ago

branches/UKMO/dev_merge_2017_restart_datestamp_GO6_mixing : recommit science changes.

File size: 44.9 KB
Line 
1MODULE zdftke
2   !!======================================================================
3   !!                       ***  MODULE  zdftke  ***
4   !! Ocean physics:  vertical mixing coefficient computed from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !! History :  OPA  !  1991-03  (b. blanke)  Original code
8   !!            7.0  !  1991-11  (G. Madec)   bug fix
9   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
10   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
11   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
12   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
13   !!            8.0  !  1997-07  (G. Madec)   lbc
14   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
15   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
16   !!             -   !  2004-10  (C. Ethe )  1D configuration
17   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
18   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
19   !!                 !           - tke penetration (wind steering)
20   !!                 !           - suface condition for tke & mixing length
21   !!                 !           - Langmuir cells
22   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
23   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
24   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
25   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
26   !!                 !                                + cleaning of the parameters + bugs correction
27   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
29   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
30   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition (ln_drg)
31   !!----------------------------------------------------------------------
32
33   !!----------------------------------------------------------------------
34   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
35   !!   tke_tke       : tke time stepping: update tke at now time step (en)
36   !!   tke_avn       : compute mixing length scale and deduce avm and avt
37   !!   zdf_tke_init  : initialization, namelist read, and parameters control
38   !!   tke_rst       : read/write tke restart in ocean restart file
39   !!----------------------------------------------------------------------
40   USE oce            ! ocean: dynamics and active tracers variables
41   USE phycst         ! physical constants
42   USE dom_oce        ! domain: ocean
43   USE domvvl         ! domain: variable volume layer
44   USE sbc_oce        ! surface boundary condition: ocean
45   USE zdfdrg         ! vertical physics: top/bottom drag coef.
46   USE zdfmxl         ! vertical physics: mixed layer
47   !
48   USE in_out_manager ! I/O manager
49   USE iom            ! I/O manager library
50   USE lib_mpp        ! MPP library
51   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
52   USE prtctl         ! Print control
53   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC   zdf_tke        ! routine called in step module
59   PUBLIC   zdf_tke_init   ! routine called in opa module
60   PUBLIC   tke_rst        ! routine called in step module
61
62   !                      !!** Namelist  namzdf_tke  **
63   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
64   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
65   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
66   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
67   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
68   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
69   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
70   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
71   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
72   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
73   LOGICAL  ::   ln_drg    ! top/bottom friction forcing flag
74   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
75   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1)
76   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
77   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
78   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells
79
80   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
81   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
82   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
83   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
84
85   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau    ! depth of tke penetration (nn_htau)
86   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl   ! now mixing lenght of dissipation
87   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr   ! now mixing lenght of dissipation
88
89   !! * Substitutions
90#  include "vectopt_loop_substitute.h90"
91   !!----------------------------------------------------------------------
92   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
93   !! $Id$
94   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
95   !!----------------------------------------------------------------------
96CONTAINS
97
98   INTEGER FUNCTION zdf_tke_alloc()
99      !!----------------------------------------------------------------------
100      !!                ***  FUNCTION zdf_tke_alloc  ***
101      !!----------------------------------------------------------------------
102      ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) ,   STAT= zdf_tke_alloc )
103      !
104      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
105      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
106      !
107   END FUNCTION zdf_tke_alloc
108
109
110   SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt )
111      !!----------------------------------------------------------------------
112      !!                   ***  ROUTINE zdf_tke  ***
113      !!
114      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
115      !!              coefficients using a turbulent closure scheme (TKE).
116      !!
117      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
118      !!              is computed from a prognostic equation :
119      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
120      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
121      !!                  + avt N^2                      ! stratif. destruc.
122      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
123      !!      with the boundary conditions:
124      !!         surface: en = max( rn_emin0, rn_ebb * taum )
125      !!         bottom : en = rn_emin
126      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
127      !!
128      !!        The now Turbulent kinetic energy is computed using the following
129      !!      time stepping: implicit for vertical diffusion term, linearized semi
130      !!      implicit for kolmogoroff dissipation term, and explicit forward for
131      !!      both buoyancy and shear production terms. Therefore a tridiagonal
132      !!      linear system is solved. Note that buoyancy and shear terms are
133      !!      discretized in a energy conserving form (Bruchard 2002).
134      !!
135      !!        The dissipative and mixing length scale are computed from en and
136      !!      the stratification (see tke_avn)
137      !!
138      !!        The now vertical eddy vicosity and diffusivity coefficients are
139      !!      given by:
140      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
141      !!              avt = max( avmb, pdl * avm                 ) 
142      !!              eav = max( avmb, avm )
143      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
144      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
145      !!
146      !! ** Action  :   compute en (now turbulent kinetic energy)
147      !!                update avt, avm (before vertical eddy coef.)
148      !!
149      !! References : Gaspar et al., JGR, 1990,
150      !!              Blanke and Delecluse, JPO, 1991
151      !!              Mellor and Blumberg, JPO 2004
152      !!              Axell, JGR, 2002
153      !!              Bruchard OM 2002
154      !!----------------------------------------------------------------------
155      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step
156      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
157      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
158      !!----------------------------------------------------------------------
159      !
160      CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt )   ! now tke (en)
161      !
162      CALL tke_avn( gdepw_n, e3t_n, e3w_n,        p_avm, p_avt )   ! now avt, avm, dissl
163      !
164  END SUBROUTINE zdf_tke
165
166
167   SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt )
168      !!----------------------------------------------------------------------
169      !!                   ***  ROUTINE tke_tke  ***
170      !!
171      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
172      !!
173      !! ** Method  : - TKE surface boundary condition
174      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
175      !!              - source term due to shear (= Kz dz[Ub] * dz[Un] )
176      !!              - Now TKE : resolution of the TKE equation by inverting
177      !!                a tridiagonal linear system by a "methode de chasse"
178      !!              - increase TKE due to surface and internal wave breaking
179      !!             NB: when sea-ice is present, both LC parameterization
180      !!                 and TKE penetration are turned off when the ice fraction
181      !!                 is smaller than 0.25
182      !!
183      !! ** Action  : - en : now turbulent kinetic energy)
184      !! ---------------------------------------------------------------------
185      USE zdf_oce , ONLY : en   ! ocean vertical physics
186      !!
187      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdepw          ! depth of w-points
188      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points)
189      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_sh2          ! shear production term
190      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
191      !
192      INTEGER ::   ji, jj, jk              ! dummy loop arguments
193      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars
194      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3
195      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient
196      REAL(wp) ::   zbbrau, zri                ! local scalars
197      REAL(wp) ::   zfact1, zfact2, zfact3     !   -         -
198      REAL(wp) ::   ztx2  , zty2  , zcof       !   -         -
199      REAL(wp) ::   ztau  , zdif               !   -         -
200      REAL(wp) ::   zus   , zwlc  , zind       !   -         -
201      REAL(wp) ::   zzd_up, zzd_lw             !   -         -
202      INTEGER , DIMENSION(jpi,jpj)     ::   imlc
203      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc
204      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
205      !!--------------------------------------------------------------------
206      !
207      zbbrau = rn_ebb / rau0       ! Local constant initialisation
208      zfact1 = -.5_wp * rdt 
209      zfact2 = 1.5_wp * rdt * rn_ediss
210      zfact3 = 0.5_wp       * rn_ediss
211      !
212      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
213      !                     !  Surface/top/bottom boundary condition on tke
214      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
215     
216      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
217         DO ji = fs_2, fs_jpim1   ! vector opt.
218            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
219         END DO
220      END DO
221      IF ( ln_isfcav ) THEN
222         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
223            DO ji = fs_2, fs_jpim1   ! vector opt.
224               en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
225            END DO
226         END DO
227      ENDIF
228      !
229      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
230      !                     !  Bottom boundary condition on tke
231      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
232      !
233      !   en(bot)   = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
234      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
235      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
236      !
237      IF( ln_drg ) THEN       !== friction used as top/bottom boundary condition on TKE
238         !
239         DO jj = 2, jpjm1           ! bottom friction
240            DO ji = fs_2, fs_jpim1     ! vector opt.
241               zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
242               zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
243               !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
244               zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
245                  &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
246               en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
247            END DO
248         END DO
249         IF( ln_isfcav ) THEN       ! top friction
250            DO jj = 2, jpjm1
251               DO ji = fs_2, fs_jpim1   ! vector opt.
252                  zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
253                  zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
254                  !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
255                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
256                     &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
257                  en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface
258               END DO
259            END DO
260         ENDIF
261         !
262      ENDIF
263      !
264      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
265      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002)
266         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
267         !
268         !                        !* total energy produce by LC : cumulative sum over jk
269         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1)
270         DO jk = 2, jpk
271            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk)
272         END DO
273         !                        !* finite Langmuir Circulation depth
274         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
275         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
276         DO jk = jpkm1, 2, -1
277            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
278               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
279                  zus  = zcof * taum(ji,jj)
280                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
281               END DO
282            END DO
283         END DO
284         !                               ! finite LC depth
285         DO jj = 1, jpj 
286            DO ji = 1, jpi
287               zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj))
288            END DO
289         END DO
290         zcof = 0.016 / SQRT( zrhoa * zcdrag )
291         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
292            DO jj = 2, jpjm1
293               DO ji = fs_2, fs_jpim1   ! vector opt.
294                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
295                  !                                           ! vertical velocity due to LC
296                  zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) )
297                  zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )
298                  !                                           ! TKE Langmuir circulation source term
299                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
300                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk)
301!!gm                     &                          / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
302               END DO
303            END DO
304         END DO
305         !
306      ENDIF
307      !
308      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
309      !                     !  Now Turbulent kinetic energy (output in en)
310      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
311      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
312      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
313      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
314      !
315      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
316         DO jk = 2, jpkm1
317            DO jj = 2, jpjm1
318               DO ji = 2, jpim1
319                  !                             ! local Richardson number
320                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
321                  !                             ! inverse of Prandtl number
322                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
323               END DO
324            END DO
325         END DO
326      ENDIF
327      !         
328      DO jk = 2, jpkm1           !* Matrix and right hand side in en
329         DO jj = 2, jpjm1
330            DO ji = fs_2, fs_jpim1   ! vector opt.
331               zcof   = zfact1 * tmask(ji,jj,jk)
332               !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
333               !                                   ! eddy coefficient (ensure numerical stability)
334               zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
335                  &          /    (  p_e3t(ji,jj,jk  ) * p_e3w(ji,jj,jk  )  )
336               zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
337                  &          /    (  p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk  )  )
338               !
339               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
340               zd_lw(ji,jj,jk) = zzd_lw
341               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
342               !
343               !                                   ! right hand side in en
344               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear
345                  &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification
346                  &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
347                  &                                ) * wmask(ji,jj,jk)
348            END DO
349         END DO
350      END DO
351      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
352      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
353         DO jj = 2, jpjm1
354            DO ji = fs_2, fs_jpim1    ! vector opt.
355               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
356            END DO
357         END DO
358      END DO
359      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
360         DO ji = fs_2, fs_jpim1   ! vector opt.
361            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
362         END DO
363      END DO
364      DO jk = 3, jpkm1
365         DO jj = 2, jpjm1
366            DO ji = fs_2, fs_jpim1    ! vector opt.
367               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)
368            END DO
369         END DO
370      END DO
371      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
372         DO ji = fs_2, fs_jpim1   ! vector opt.
373            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
374         END DO
375      END DO
376      DO jk = jpk-2, 2, -1
377         DO jj = 2, jpjm1
378            DO ji = fs_2, fs_jpim1    ! vector opt.
379               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
380            END DO
381         END DO
382      END DO
383      DO jk = 2, jpkm1                             ! set the minimum value of tke
384         DO jj = 2, jpjm1
385            DO ji = fs_2, fs_jpim1   ! vector opt.
386               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
387            END DO
388         END DO
389      END DO
390      !
391      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
392      !                            !  TKE due to surface and internal wave breaking
393      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
394!!gm BUG : in the exp  remove the depth of ssh !!!
395!!gm       i.e. use gde3w in argument (pdepw)
396     
397     
398      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
399         DO jk = 2, jpkm1
400            DO jj = 2, jpjm1
401               DO ji = fs_2, fs_jpim1   ! vector opt.
402                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
403                     &                                 * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk)
404!!gm                     &                             * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
405!!gm    multiplication by surface tmask useless (already includes in MAX( 0, 1-4*fr_i )
406               END DO
407            END DO
408         END DO
409      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
410         DO jj = 2, jpjm1
411            DO ji = fs_2, fs_jpim1   ! vector opt.
412               jk = nmln(ji,jj)
413               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
414                  &                                 * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk)
415!!gm                  &                             * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
416            END DO
417         END DO
418      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
419         DO jk = 2, jpkm1
420            DO jj = 2, jpjm1
421               DO ji = fs_2, fs_jpim1   ! vector opt.
422                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
423                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
424                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
425                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
426                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
427                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
428                     &                        * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk)
429!!gm                     &                    * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
430               END DO
431            END DO
432         END DO
433      ENDIF
434      !
435   END SUBROUTINE tke_tke
436
437
438   SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt )
439      !!----------------------------------------------------------------------
440      !!                   ***  ROUTINE tke_avn  ***
441      !!
442      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
443      !!
444      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
445      !!              the tke_tke routine). First, the now mixing lenth is
446      !!      computed from en and the strafification (N^2), then the mixings
447      !!      coefficients are computed.
448      !!              - Mixing length : a first evaluation of the mixing lengh
449      !!      scales is:
450      !!                      mxl = sqrt(2*en) / N 
451      !!      where N is the brunt-vaisala frequency, with a minimum value set
452      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
453      !!        The mixing and dissipative length scale are bound as follow :
454      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
455      !!                        zmxld = zmxlm = mxl
456      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
457      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
458      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
459      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
460      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
461      !!                    the surface to obtain ldown. the resulting length
462      !!                    scales are:
463      !!                        zmxld = sqrt( lup * ldown )
464      !!                        zmxlm = min ( lup , ldown )
465      !!              - Vertical eddy viscosity and diffusivity:
466      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
467      !!                      avt = max( avmb, pdlr * avm ) 
468      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
469      !!
470      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
471      !!----------------------------------------------------------------------
472      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics
473      !!
474      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! depth (w-points)
475      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points)
476      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
477      !
478      INTEGER  ::   ji, jj, jk   ! dummy loop indices
479      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
480      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
481      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      -
482      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace
483      !!--------------------------------------------------------------------
484      !
485      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
486      !                     !  Mixing length
487      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
488      !
489      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
490      !
491      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
492      zmxlm(:,:,:)  = rmxl_min   
493      zmxld(:,:,:)  = rmxl_min
494      !
495      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
496         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
497         DO jj = 2, jpjm1
498            DO ji = fs_2, fs_jpim1
499               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
500            END DO
501         END DO
502      ELSE
503         zmxlm(:,:,1) = rn_mxl0
504      ENDIF
505      !
506      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
507         DO jj = 2, jpjm1
508            DO ji = fs_2, fs_jpim1   ! vector opt.
509               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
510               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
511            END DO
512         END DO
513      END DO
514      !
515      !                     !* Physical limits for the mixing length
516      !
517      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
518      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
519      !
520      SELECT CASE ( nn_mxl )
521      !
522 !!gm Not sure of that coding for ISF....
523      ! where wmask = 0 set zmxlm == p_e3w
524      CASE ( 0 )           ! bounded by the distance to surface and bottom
525         DO jk = 2, jpkm1
526            DO jj = 2, jpjm1
527               DO ji = fs_2, fs_jpim1   ! vector opt.
528                  zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
529                  &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) )
530                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
531                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
532                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
533               END DO
534            END DO
535         END DO
536         !
537      CASE ( 1 )           ! bounded by the vertical scale factor
538         DO jk = 2, jpkm1
539            DO jj = 2, jpjm1
540               DO ji = fs_2, fs_jpim1   ! vector opt.
541                  zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) )
542                  zmxlm(ji,jj,jk) = zemxl
543                  zmxld(ji,jj,jk) = zemxl
544               END DO
545            END DO
546         END DO
547         !
548      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
549         DO jk = 2, jpkm1         ! from the surface to the bottom :
550            DO jj = 2, jpjm1
551               DO ji = fs_2, fs_jpim1   ! vector opt.
552                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
553               END DO
554            END DO
555         END DO
556         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
557            DO jj = 2, jpjm1
558               DO ji = fs_2, fs_jpim1   ! vector opt.
559                  zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
560                  zmxlm(ji,jj,jk) = zemxl
561                  zmxld(ji,jj,jk) = zemxl
562               END DO
563            END DO
564         END DO
565         !
566      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
567         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
568            DO jj = 2, jpjm1
569               DO ji = fs_2, fs_jpim1   ! vector opt.
570                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
571               END DO
572            END DO
573         END DO
574         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
575            DO jj = 2, jpjm1
576               DO ji = fs_2, fs_jpim1   ! vector opt.
577                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
578               END DO
579            END DO
580         END DO
581         DO jk = 2, jpkm1
582            DO jj = 2, jpjm1
583               DO ji = fs_2, fs_jpim1   ! vector opt.
584                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
585                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
586                  zmxlm(ji,jj,jk) = zemlm
587                  zmxld(ji,jj,jk) = zemlp
588               END DO
589            END DO
590         END DO
591         !
592      END SELECT
593      !
594      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
595      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
596      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
597      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
598         DO jj = 2, jpjm1
599            DO ji = fs_2, fs_jpim1   ! vector opt.
600               zsqen = SQRT( en(ji,jj,jk) )
601               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
602               p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
603               p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
604               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
605            END DO
606         END DO
607      END DO
608      !
609      !
610      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
611         DO jk = 2, jpkm1
612            DO jj = 2, jpjm1
613               DO ji = fs_2, fs_jpim1   ! vector opt.
614                  p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
615              END DO
616            END DO
617         END DO
618      ENDIF
619      !
620      IF(ln_ctl) THEN
621         CALL prt_ctl( tab3d_1=en   , clinfo1=' tke  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
622         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke  - m: ', kdim=jpk )
623      ENDIF
624      !
625   END SUBROUTINE tke_avn
626
627
628   SUBROUTINE zdf_tke_init
629      !!----------------------------------------------------------------------
630      !!                  ***  ROUTINE zdf_tke_init  ***
631      !!                     
632      !! ** Purpose :   Initialization of the vertical eddy diffivity and
633      !!              viscosity when using a tke turbulent closure scheme
634      !!
635      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
636      !!              called at the first timestep (nit000)
637      !!
638      !! ** input   :   Namlist namzdf_tke
639      !!
640      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
641      !!----------------------------------------------------------------------
642      USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag
643      !!
644      INTEGER ::   ji, jj, jk   ! dummy loop indices
645      INTEGER ::   ios
646      !!
647      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,          &
648         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,          &
649         &                 rn_mxl0 , nn_pdl   , ln_drg , ln_lc    , rn_lc,   &
650         &                 nn_etau , nn_htau  , rn_efr   
651      !!----------------------------------------------------------------------
652      !
653      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
654      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
655901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
656
657      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
658      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
659902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
660      IF(lwm) WRITE ( numond, namzdf_tke )
661      !
662      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
663      !
664      IF(lwp) THEN                    !* Control print
665         WRITE(numout,*)
666         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
667         WRITE(numout,*) '~~~~~~~~~~~~'
668         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
669         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
670         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
671         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
672         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
673         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
674         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
675         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
676         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
677         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
678         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
679         WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg
680         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
681         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
682         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
683         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
684         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr
685         IF( ln_drg ) THEN
686            WRITE(numout,*)
687            WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
688            WRITE(numout,*) '      top    ocean cavity roughness (m)          rn_z0(_top)= ', r_z0_top
689            WRITE(numout,*) '      Bottom seafloor     roughness (m)          rn_z0(_bot)= ', r_z0_bot
690         ENDIF
691         WRITE(numout,*)
692         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
693         WRITE(numout,*)
694      ENDIF
695      !
696      IF( ln_zdfiwm ) THEN          ! Internal wave-driven mixing
697         rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
698         rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
699         IF(lwp) WRITE(numout,*) '   ==>>>   Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3'
700      ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
701         rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
702         IF(lwp) WRITE(numout,*) '   ==>>>   minimum mixing length with your parameters rmxl_min = ', rmxl_min
703      ENDIF
704      !
705      !                              ! allocate tke arrays
706      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
707      !
708      !                               !* Check of some namelist values
709      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
710      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
711      IF( nn_htau < 0 .OR. nn_htau == 2 .OR. nn_htau > 5 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1, 3, 4 or 5 ' )
712      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
713      !
714      IF( ln_mxl0 ) THEN
715         IF(lwp) WRITE(numout,*)
716         IF(lwp) WRITE(numout,*) '   ==>>>   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
717         rn_mxl0 = rmxl_min
718      ENDIF
719     
720      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
721
722      !                               !* depth of penetration of surface tke
723      IF( nn_etau /= 0 ) THEN     
724         htau(:,:) = 0._wp
725         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
726         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
727            htau(:,:) = 10._wp
728         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
729            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
730         CASE( 3 )                                 ! F(latitude) : 0.5m to 15m poleward of 20 degrees
731            htau(:,:) = MAX(  0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
732         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south
733            DO jj = 2, jpjm1
734               DO ji = fs_2, fs_jpim1   ! vector opt.
735                  IF( gphit(ji,jj) <= 0._wp ) THEN
736                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
737                  ELSE
738                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
739                  ENDIF
740               END DO
741            END DO
742         CASE ( 5 )                                ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south,
743            DO jj = 2, jpjm1                       !               10m to 30m between 30/45 degrees south
744               DO ji = fs_2, fs_jpim1   ! vector opt.
745                  IF( gphit(ji,jj) <= -30._wp ) THEN
746                     htau(ji,jj) = MAX(  10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) )   )
747                  ELSE
748                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
749                  ENDIF
750               END DO
751            END DO
752         END SELECT
753         !
754      ENDIF
755      !                                !* read or initialize all required files
756      CALL tke_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, dissl)
757      !
758      IF( lwxios ) THEN
759         CALL iom_set_rstw_var_active('en')
760         CALL iom_set_rstw_var_active('avt_k')
761         CALL iom_set_rstw_var_active('avm_k')
762         CALL iom_set_rstw_var_active('dissl')
763      ENDIF
764   END SUBROUTINE zdf_tke_init
765
766
767   SUBROUTINE tke_rst( kt, cdrw )
768      !!---------------------------------------------------------------------
769      !!                   ***  ROUTINE tke_rst  ***
770      !!                     
771      !! ** Purpose :   Read or write TKE file (en) in restart file
772      !!
773      !! ** Method  :   use of IOM library
774      !!                if the restart does not contain TKE, en is either
775      !!                set to rn_emin or recomputed
776      !!----------------------------------------------------------------------
777      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
778      !!
779      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
780      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
781      !
782      INTEGER ::   jit, jk              ! dummy loop indices
783      INTEGER ::   id1, id2, id3, id4   ! local integers
784      !!----------------------------------------------------------------------
785      !
786      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
787         !                                   ! ---------------
788         IF( ln_rstart ) THEN                   !* Read the restart file
789            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
790            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
791            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
792            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
793            !
794            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
795               CALL iom_get( numror, jpdom_autoglo, 'en'   , en   , ldxios = lrxios )
796               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k, ldxios = lrxios )
797               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k, ldxios = lrxios )
798               CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl, ldxios = lrxios )
799            ELSE                                          ! start TKE from rest
800               IF(lwp) WRITE(numout,*)
801               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without TKE scheme, set en to background values'
802               en   (:,:,:) = rn_emin * wmask(:,:,:)
803               dissl(:,:,:) = 1.e-12_wp
804               ! avt_k, avm_k already set to the background value in zdf_phy_init
805            ENDIF
806         ELSE                                   !* Start from rest
807            IF(lwp) WRITE(numout,*)
808            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set en to the background value'
809            en   (:,:,:) = rn_emin * wmask(:,:,:)
810            dissl(:,:,:) = 1.e-12_wp
811            ! avt_k, avm_k already set to the background value in zdf_phy_init
812         ENDIF
813         !
814      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
815         !                                   ! -------------------
816         IF(lwp) WRITE(numout,*) '---- tke_rst ----'
817         IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
818         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en   , ldxios = lwxios )
819         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios )
820         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios )
821         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl, ldxios = lwxios )
822         IF( lwxios ) CALL iom_swap(      cxios_context          )
823         !
824      ENDIF
825      !
826   END SUBROUTINE tke_rst
827
828   !!======================================================================
829END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.