source: NEMO/trunk/src/OCE/ZDF/zdfgls.F90 @ 14966

Last change on this file since 14966 was 14966, checked in by cbricaud, 8 months ago

add scaling of mixing under sea-ice in GLS (as for TKE) ; see ticket #2604

  • Property svn:keywords set to Id
File size: 61.0 KB
Line 
1MODULE zdfgls
2   !!======================================================================
3   !!                       ***  MODULE  zdfgls  ***
4   !! Ocean physics:  vertical mixing coefficient computed from the gls
5   !!                 turbulent closure parameterization
6   !!======================================================================
7   !! History :  3.0  !  2009-09  (G. Reffray)  Original code
8   !!            3.3  !  2010-10  (C. Bricaud)  Add in the reference
9   !!            4.0  !  2017-04  (G. Madec)  remove CPP keys & avm at t-point only
10   !!             -   !  2017-05  (G. Madec)  add top friction as boundary condition
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   zdf_gls       : update momentum and tracer Kz from a gls scheme
15   !!   zdf_gls_init  : initialization, namelist read, and parameters control
16   !!   gls_rst       : read/write gls restart in ocean restart file
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and active tracers
19   USE dom_oce        ! ocean space and time domain
20   USE domvvl         ! ocean space and time domain : variable volume layer
21   USE zdfdrg  , ONLY : ln_drg_OFF            ! top/bottom free-slip flag
22   USE zdfdrg  , ONLY : r_z0_top , r_z0_bot   ! top/bottom roughness
23   USE zdfdrg  , ONLY : rCdU_top , rCdU_bot   ! top/bottom friction
24   USE sbc_oce        ! surface boundary condition: ocean
25   USE phycst         ! physical constants
26   USE zdfmxl         ! mixed layer
27   USE sbcwave , ONLY : hsw   ! significant wave height
28#if defined key_si3
29   USE ice, ONLY: hm_i, h_i
30#endif
31#if defined key_cice
32   USE sbc_ice, ONLY: h_i
33#endif
34   !
35   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
36   USE lib_mpp        ! MPP manager
37   USE prtctl         ! Print control
38   USE in_out_manager ! I/O manager
39   USE iom            ! I/O manager library
40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   zdf_gls        ! called in zdfphy
46   PUBLIC   zdf_gls_init   ! called in zdfphy
47   PUBLIC   gls_rst        ! called in zdfphy
48
49   !
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxl_n    !: now mixing length
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_surf !: Squared surface velocity scale at T-points
53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_top  !: Squared top     velocity scale at T-points
54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_bot  !: Squared bottom  velocity scale at T-points
55
56   !                              !! ** Namelist  namzdf_gls  **
57   LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
58   LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
59   INTEGER  ::   nn_mxlice         ! type of scaling under sea-ice (=0/1/2/3)
60   INTEGER  ::   nn_bc_surf        ! surface boundary condition (=0/1)
61   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1)
62   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation
63   INTEGER  ::   nn_z0_ice         ! Roughness accounting for sea ice
64   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
65   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
66   REAL(wp) ::   rn_clim_galp      ! Holt 2008 value for k-eps: 0.267
67   REAL(wp) ::   rn_epsmin         ! minimum value of dissipation (m2/s3)
68   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
69   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
70   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
71   REAL(wp) ::   rn_hsro           ! Minimum surface roughness
72   REAL(wp) ::   rn_hsri           ! Ice ocean roughness
73   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)
74
75   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
76   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
77   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn
78   REAL(wp) ::   rghmin        = -0.28_wp
79   REAL(wp) ::   rgh0          =  0.0329_wp
80   REAL(wp) ::   rghcri        =  0.03_wp
81   REAL(wp) ::   ra1           =  0.92_wp
82   REAL(wp) ::   ra2           =  0.74_wp
83   REAL(wp) ::   rb1           = 16.60_wp
84   REAL(wp) ::   rb2           = 10.10_wp
85   REAL(wp) ::   re2           =  1.33_wp
86   REAL(wp) ::   rl1           =  0.107_wp
87   REAL(wp) ::   rl2           =  0.0032_wp
88   REAL(wp) ::   rl3           =  0.0864_wp
89   REAL(wp) ::   rl4           =  0.12_wp
90   REAL(wp) ::   rl5           = 11.9_wp
91   REAL(wp) ::   rl6           =  0.4_wp
92   REAL(wp) ::   rl7           =  0.0_wp
93   REAL(wp) ::   rl8           =  0.48_wp
94   REAL(wp) ::   rm1           =  0.127_wp
95   REAL(wp) ::   rm2           =  0.00336_wp
96   REAL(wp) ::   rm3           =  0.0906_wp
97   REAL(wp) ::   rm4           =  0.101_wp
98   REAL(wp) ::   rm5           = 11.2_wp
99   REAL(wp) ::   rm6           =  0.4_wp
100   REAL(wp) ::   rm7           =  0.0_wp
101   REAL(wp) ::   rm8           =  0.318_wp
102   REAL(wp) ::   rtrans        =  0.1_wp
103   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
104   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        -
105   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        -
106   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
107   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
108   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        -
109   REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5                     !     -           -           -        -
110   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        -
111   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        -
112   !
113   REAL(wp) ::   r2_3 = 2._wp/3._wp   ! constant=2/3
114
115   !! * Substitutions
116#  include "do_loop_substitute.h90"
117#  include "domzgr_substitute.h90"
118   !!----------------------------------------------------------------------
119   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
120   !! $Id$
121   !! Software governed by the CeCILL license (see ./LICENSE)
122   !!----------------------------------------------------------------------
123CONTAINS
124
125   INTEGER FUNCTION zdf_gls_alloc()
126      !!----------------------------------------------------------------------
127      !!                ***  FUNCTION zdf_gls_alloc  ***
128      !!----------------------------------------------------------------------
129      ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) ,                     &
130         &      zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc )
131         !
132      CALL mpp_sum ( 'zdfgls', zdf_gls_alloc )
133      IF( zdf_gls_alloc /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_alloc: failed to allocate arrays' )
134   END FUNCTION zdf_gls_alloc
135
136
137   SUBROUTINE zdf_gls( kt, Kbb, Kmm, p_sh2, p_avm, p_avt )
138      !!----------------------------------------------------------------------
139      !!                   ***  ROUTINE zdf_gls  ***
140      !!
141      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
142      !!              coefficients using the GLS turbulent closure scheme.
143      !!----------------------------------------------------------------------
144      USE zdf_oce , ONLY : en, avtb, avmb   ! ocean vertical physics
145      !!
146      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time step
147      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
148      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
149      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
150      !
151      INTEGER  ::   ji, jj, jk    ! dummy loop arguments
152      INTEGER  ::   ibot, ibotm1  ! local integers
153      INTEGER  ::   itop, itopp1  !   -       -
154      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1 , zex2  ! local scalars
155      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof, zdir  !   -      -
156      REAL(wp) ::   zratio, zrn2, zflxb, sh     , z_en  !   -      -
157      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
158      REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      -
159      REAL(wp) ::   zmsku, zmskv                        !   -      -
160      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zdep
161      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zkar
162      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zflxs                 ! Turbulence fluxed induced by internal waves
163      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zhsro                 ! Surface roughness (surface waves)
164      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra              ! Tapering of wave breaking under sea ice
165      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eb                    ! tke at time before
166      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   hmxl_b                ! mixing length at time before
167      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eps                   ! dissipation rate
168      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwall_psi             ! Wall function use in the wb case (ln_sigpsi)
169      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   psi                   ! psi at time now
170      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
171      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zstt, zstm            ! stability function on tracer and momentum
172      !!--------------------------------------------------------------------
173      !
174      ! Preliminary computing
175      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
176         ustar2_surf(ji,jj) = 0._wp   ;   ustar2_top(ji,jj) = 0._wp   ;   ustar2_bot(ji,jj) = 0._wp
177      END_2D
178      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
179         psi(ji,jj,jk) = 0._wp   ;   zwall_psi(ji,jj,jk) = 0._wp
180      END_3D
181
182      SELECT CASE ( nn_z0_ice )
183      CASE( 0 )   ;   zice_fra(:,:) = 0._wp
184      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(A2D(nn_hls)) * 10._wp )
185      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(A2D(nn_hls))
186      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
187      END SELECT
188
189      ! Compute surface, top and bottom friction at T-points
190      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          !==  surface ocean friction
191         ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction
192      END_2D
193      !
194      !!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that...
195      !
196      IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction)
197         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          ! bottom friction (explicit before friction)
198            zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
199            zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
200            ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  &
201               &                                         + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
202         END_2D
203         IF( ln_isfcav ) THEN
204            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )      ! top friction
205               zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
206               zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
207               ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  &
208                  &                                         + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
209            END_2D
210         ENDIF
211      ENDIF
212
213      SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==!
214      CASE ( 0 )                          ! Constant roughness
215         zhsro(:,:) = rn_hsro
216      CASE ( 1 )             ! Standard Charnock formula
217         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
218            zhsro(ji,jj) = MAX( rsbc_zs1 * ustar2_surf(ji,jj) , rn_hsro )
219         END_2D
220      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
221!!gm faster coding : the 2 comment lines should be used
222!!gm         zcof = 2._wp * 0.6_wp / 28._wp
223!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) )  )       ! Wave age (eq. 10)
224         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
225            zcof = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(ji,jj),rsmall))) )          ! Wave age (eq. 10)
226            zhsro(ji,jj) = MAX(rsbc_zs2 * ustar2_surf(ji,jj) * zcof**1.5, rn_hsro)        ! zhsro = rn_frac_hs * Hsw (eq. 11)
227         END_2D
228      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
229         zhsro(:,:) = MAX(rn_frac_hs * hsw(A2D(nn_hls)), rn_hsro)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )
230      END SELECT
231      !
232      ! adapt roughness where there is sea ice
233      SELECT CASE( nn_mxlice )       ! Type of scaling under sea-ice
234      !
235      CASE( 1 )                      ! scaling with constant sea-ice roughness
236         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
237            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro
238         END_2D
239         !
240      CASE( 2 )                      ! scaling with mean sea-ice thickness
241#if defined key_si3
242         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
243            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * hm_i(ji,jj) )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro
244         END_2D
245#endif
246         !
247      CASE( 3 )                      ! scaling with max sea-ice thickness
248#if defined key_si3 || defined key_cice
249         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
250            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * MAXVAL(h_i(ji,jj,:)) )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro
251         END_2D
252#endif
253         !
254      END SELECT
255      !
256      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  !==  Compute dissipation rate  ==!
257         eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
258      END_3D
259
260      ! Save tke at before time step
261      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
262         eb    (ji,jj,jk) = en    (ji,jj,jk)
263         hmxl_b(ji,jj,jk) = hmxl_n(ji,jj,jk)
264      END_3D
265
266      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
267         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
268            zup   = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
269            zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )
270            zcoef = ( zup / MAX( zdown, rsmall ) )
271            zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk)
272         END_3D
273      ENDIF
274
275      !!---------------------------------!!
276      !!   Equation to prognostic k      !!
277      !!---------------------------------!!
278      !
279      ! Now Turbulent kinetic energy (output in en)
280      ! -------------------------------
281      ! Resolution of a tridiagonal linear system by a "methode de chasse"
282      ! computation from level 2 to jpkm1  (e(1) computed after and e(jpk)=0 ).
283      ! The surface boundary condition are set after
284      ! The bottom boundary condition are also set after. In standard e(bottom)=0.
285      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
286      ! Warning : after this step, en : right hand side of the matrix
287
288      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
289         !
290         buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
291         !
292         diss = eps(ji,jj,jk)                         ! dissipation
293         !
294         zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy )   ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
295         !
296         zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk)          ! production term
297         zdiss = zdir*(diss/en(ji,jj,jk))   +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
298!!gm better coding, identical results
299!               zesh2 =   p_sh2(ji,jj,jk) + zdir*buoy               ! production term
300!               zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term
301!!gm
302         !
303         ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
304         ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style)
305         ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
306         ! Otherwise, this should be rsc_psi/rsc_psi0
307         IF( ln_sigpsi ) THEN
308            zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
309            zwall_psi(ji,jj,jk) = rsc_psi /   &
310               &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
311         ELSE
312            zwall_psi(ji,jj,jk) = 1._wp
313         ENDIF
314         !
315         ! building the matrix
316         zcof = rfact_tke * tmask(ji,jj,jk)
317         !                                        ! lower diagonal, in fact not used for jk = 2 (see surface conditions)
318         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) )   &
319            &                 / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
320         !                                        ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions)
321         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) )   &
322            &                 / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
323         !                                        ! diagonal
324         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rn_Dt * zdiss * wmask(ji,jj,jk)
325         !                                        ! right hand side in en
326         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk)
327      END_3D
328      !
329      zdiag(:,:,jpk) = 1._wp
330      !
331      ! Set surface condition on zwall_psi (1 at the bottom)
332      zwall_psi(:,:, 1 ) = zwall_psi(:,:,2)
333      zwall_psi(:,:,jpk) = 1._wp
334      !
335      ! Surface boundary condition on tke
336      ! ---------------------------------
337      !
338      SELECT CASE ( nn_bc_surf )
339      !
340      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2)
341         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
342            ! First level
343            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  )
344            zd_lw(ji,jj,1) = en(ji,jj,1)
345            zd_up(ji,jj,1) = 0._wp
346            zdiag(ji,jj,1) = 1._wp
347            !
348            ! One level below
349            en   (ji,jj,2) =  MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1          &
350               &                             * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )**(1.5_wp*ra_sf)  )**r2_3 )
351            zd_lw(ji,jj,2) = 0._wp
352            zd_up(ji,jj,2) = 0._wp
353            zdiag(ji,jj,2) = 1._wp
354         END_2D
355         !
356         !
357      CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
358         !
359         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
360            ! Dirichlet conditions at k=1
361            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  )
362            zd_lw(ji,jj,1) = en(ji,jj,1)
363            zd_up(ji,jj,1) = 0._wp
364            zdiag(ji,jj,1) = 1._wp
365            !
366            ! at k=2, set de/dz=Fw
367            !cbr
368            ! zdiag zd_lw not defined/used on the halo
369            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
370            zd_lw(ji,jj,2) = 0._wp
371            !
372            zkar (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) ))
373            zflxs(ji,jj)   = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) &
374                &                    * (  ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj)  )**(1.5_wp*ra_sf)
375!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf)
376            en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
377         END_2D
378         !
379         !
380      END SELECT
381
382      ! Bottom boundary condition on tke
383      ! --------------------------------
384      !
385      SELECT CASE ( nn_bc_bot )
386      !
387      CASE ( 0 )             ! Dirichlet
388         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
389         !                      ! Balance between the production and the dissipation terms
390         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
391!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
392!!   With thick deep ocean level thickness, this may be quite large, no ???
393!!   in particular in ocean cavities where top stratification can be large...
394            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
395            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
396            !
397            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
398            !
399            ! Dirichlet condition applied at:
400            !     Bottom level (ibot)      &      Just above it (ibotm1)
401            zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp
402            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
403            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp
404            en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en
405         END_2D
406         !
407         ! NOTE: ctl_stop with ln_isfcav when using GLS
408         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
409            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
410               itop   = mikt(ji,jj)       ! k   top w-point
411               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
412               !                                                ! mask at the ocean surface points
413               z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
414               !
415 !!gm TO BE VERIFIED !!!
416               ! Dirichlet condition applied at:
417               !     top level (itop)         &      Just below it (itopp1)
418               zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
419               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
420               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
421               en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
422            END_2D
423         ENDIF
424         !
425      CASE ( 1 )             ! Neumman boundary condition
426         !
427         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
428            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
429            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
430            !
431            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
432            !
433            ! Bottom level Dirichlet condition:
434            !     Bottom level (ibot)      &      Just above it (ibotm1)
435            !         Dirichlet            !         Neumann
436            zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag
437            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
438            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
439            en   (ji,jj,ibot) = z_en
440         END_2D
441         ! NOTE: ctl_stop with ln_isfcav when using GLS
442         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
443            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
444               itop   = mikt(ji,jj)       ! k   top w-point
445               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
446               !                                                ! mask at the ocean surface points
447               z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
448               !
449               ! Bottom level Dirichlet condition:
450               !     Bottom level (ibot)      &      Just above it (ibotm1)
451               !         Dirichlet            !         Neumann
452               zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
453               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
454               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
455               en   (ji,jj,itop) = z_en
456            END_2D
457         ENDIF
458         !
459      END SELECT
460
461      ! Matrix inversion (en prescribed at surface and the bottom)
462      ! ----------------------------------------------------------
463      !
464      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
465         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
466      END_3D
467      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
468         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)
469      END_3D
470      DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
471         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
472      END_3D
473      !                                            ! set the minimum value of tke
474      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
475         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin )
476      END_3D
477
478      !!----------------------------------------!!
479      !!   Solve prognostic equation for psi    !!
480      !!----------------------------------------!!
481
482      ! Set psi to previous time step value
483      !
484      SELECT CASE ( nn_clos )
485      !
486      CASE( 0 )               ! k-kl  (Mellor-Yamada)
487         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
488            psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
489         END_3D
490         !
491      CASE( 1 )               ! k-eps
492         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
493            psi(ji,jj,jk)  = eps(ji,jj,jk)
494         END_3D
495         !
496      CASE( 2 )               ! k-w
497         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
498            psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
499         END_3D
500         !
501      CASE( 3 )               ! generic
502         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
503            psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn
504         END_3D
505         !
506      END SELECT
507      !
508      ! Now gls (output in psi)
509      ! -------------------------------
510      ! Resolution of a tridiagonal linear system by a "methode de chasse"
511      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
512      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
513      ! Warning : after this step, en : right hand side of the matrix
514
515      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
516         !
517         ! psi / k
518         zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
519         !
520         ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable)
521         zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
522         !
523         rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p
524         !
525         ! shear prod. - stratif. destruction
526         prod = rpsi1 * zratio * p_sh2(ji,jj,jk)
527         !
528         ! stratif. destruction
529         buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
530         !
531         ! shear prod. - stratif. destruction
532         diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
533         !
534         zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
535         !
536         zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term
537         zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
538         !
539         ! building the matrix
540         zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
541         !                                               ! lower diagonal
542         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) )   &
543            &            / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
544         !                                               ! upper diagonal
545         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) )   &
546            &            / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
547         !                                               ! diagonal
548         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk)
549         !                                               ! right hand side in psi
550         psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk)
551      END_3D
552      !
553      zdiag(:,:,jpk) = 1._wp
554
555      ! Surface boundary condition on psi
556      ! ---------------------------------
557      !
558      SELECT CASE ( nn_bc_surf )
559      !
560      CASE ( 0 )             ! Dirichlet boundary conditions
561         !
562         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
563            ! Surface value
564            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf ! Cosmetic
565            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
566            zd_lw(ji,jj,1) = psi(ji,jj,1)
567            zd_up(ji,jj,1) = 0._wp
568            zdiag(ji,jj,1) = 1._wp
569            !
570            ! One level below
571            zkar    (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) )))
572            zdep    (ji,jj)   = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj)
573            psi     (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
574            zd_lw(ji,jj,2) = 0._wp
575            zd_up(ji,jj,2) = 0._wp
576            zdiag(ji,jj,2) = 1._wp
577         END_2D
578         !
579      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
580         !
581         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
582            ! Surface value: Dirichlet
583            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf
584            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
585            zd_lw(ji,jj,1) = psi(ji,jj,1)
586            zd_up(ji,jj,1) = 0._wp
587            zdiag(ji,jj,1) = 1._wp
588            !
589            ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo
590            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
591            zd_lw(ji,jj,2) = 0._wp
592            !
593            ! Set psi vertical flux at the surface:
594            zkar (ji,jj)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope
595            zdep (ji,jj)   = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf)
596            zflxs(ji,jj)   = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) &
597               &           *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp)
598            zdep (ji,jj)   = rsbc_psi1 * (zwall_psi(ji,jj,1)*p_avm(ji,jj,1)+zwall_psi(ji,jj,2)*p_avm(ji,jj,2)) * &
599               &           ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.)
600            zflxs(ji,jj)   = zdep(ji,jj) * zflxs(ji,jj)
601            psi  (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
602         END_2D
603         !
604      END SELECT
605
606      ! Bottom boundary condition on psi
607      ! --------------------------------
608      !
609!!gm should be done for ISF (top boundary cond.)
610!!gm so, totally new staff needed      ===>>> think about that !
611!
612      SELECT CASE ( nn_bc_bot )     ! bottom boundary
613      !
614      CASE ( 0 )             ! Dirichlet
615         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
616         !                      ! Balance between the production and the dissipation terms
617         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
618            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
619            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
620            zdep(ji,jj) = vkarmn * r_z0_bot
621            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
622            zd_lw(ji,jj,ibot) = 0._wp
623            zd_up(ji,jj,ibot) = 0._wp
624            zdiag(ji,jj,ibot) = 1._wp
625            !
626            ! Just above last level, Dirichlet condition again (GOTM like)
627            zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) )
628            psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
629            zd_lw(ji,jj,ibotm1) = 0._wp
630            zd_up(ji,jj,ibotm1) = 0._wp
631            zdiag(ji,jj,ibotm1) = 1._wp
632         END_2D
633         !
634      CASE ( 1 )             ! Neumman boundary condition
635         !
636         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
637            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
638            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
639            !
640            ! Bottom level Dirichlet condition:
641            zdep(ji,jj) = vkarmn * r_z0_bot
642            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
643            !
644            zd_lw(ji,jj,ibot) = 0._wp
645            zd_up(ji,jj,ibot) = 0._wp
646            zdiag(ji,jj,ibot) = 1._wp
647            !
648            ! Just above last level: Neumann condition with flux injection
649            zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag
650            zd_up(ji,jj,ibotm1) = 0.
651            !
652            ! Set psi vertical flux at the bottom:
653            zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm)
654            zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
655               &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
656            psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm)
657         END_2D
658         !
659      END SELECT
660
661      ! Matrix inversion
662      ! ----------------
663      !
664      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
665         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
666      END_3D
667      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
668         zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
669      END_3D
670      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
671         psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
672      END_3D
673
674      ! Set dissipation
675      !----------------
676
677      SELECT CASE ( nn_clos )
678      !
679      CASE( 0 )               ! k-kl  (Mellor-Yamada)
680         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
681            eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
682         END_3D
683         !
684      CASE( 1 )               ! k-eps
685         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
686            eps(ji,jj,jk) = psi(ji,jj,jk)
687         END_3D
688         !
689      CASE( 2 )               ! k-w
690         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
691            eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
692         END_3D
693         !
694      CASE( 3 )               ! generic
695         zcoef = rc0**( 3._wp  + rpp/rnn )
696         zex1  =      ( 1.5_wp + rmm/rnn )
697         zex2  = -1._wp / rnn
698         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
699            eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
700         END_3D
701         !
702      END SELECT
703
704      ! Limit dissipation rate under stable stratification
705      ! --------------------------------------------------
706      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
707         ! limitation
708         eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
709         hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
710      END_3D
711      IF( ln_length_lim ) THEN        ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
712         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
713            zrn2 = MAX( rn2(ji,jj,jk), rsmall )
714            hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
715         END_3D
716      ENDIF
717
718      !
719      ! Stability function and vertical viscosity and diffusivity
720      ! ---------------------------------------------------------
721      !
722      SELECT CASE ( nn_stab_func )
723      !
724      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
725         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
726            ! zcof =  l²/q²
727            zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
728            ! Gh = -N²l²/q²
729            gh = - rn2(ji,jj,jk) * zcof
730            gh = MIN( gh, rgh0   )
731            gh = MAX( gh, rghmin )
732            ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
733            sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) )
734            sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh)
735            !
736            ! Store stability function in zstt and zstm
737            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
738            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
739         END_3D
740         !
741      CASE ( 2, 3 )               ! Canuto stability functions
742         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
743            ! zcof =  l²/q²
744            zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
745            ! Gh = -N²l²/q²
746            gh = - rn2(ji,jj,jk) * zcof
747            gh = MIN( gh, rgh0   )
748            gh = MAX( gh, rghmin )
749            gh = gh * rf6
750            ! Gm =  M²l²/q² Shear number
751            shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
752            gm = MAX( shr * zcof , 1.e-10 )
753            gm = gm * rf6
754            gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm )
755            ! Stability functions from Canuto
756            rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm
757            sm = (rs0 - rs1*gh + rs2*gm) / rcff
758            sh = (rs4 - rs5*gh + rs6*gm) / rcff
759            !
760            ! Store stability function in zstt and zstm
761            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
762            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
763         END_3D
764         !
765      END SELECT
766
767      ! Boundary conditions on stability functions for momentum (Neumann):
768      ! Lines below are useless if GOTM style Dirichlet conditions are used
769
770      zstm(:,:,1) = zstm(:,:,2)
771
772      ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
773      zstm(:,:,jpk) = 0.
774      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )             ! update bottom with good values
775         zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
776      END_2D
777
778      zstt(:,:,  1) = wmask(A2D(nn_hls),  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
779      zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
780
781!!gm should be done for ISF (top boundary cond.)
782!!gm so, totally new staff needed!!gm
783
784      ! Compute diffusivities/viscosities
785      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
786      !  -> yes BUT p_avm(:,:1) and p_avm(:,:jpk) are used when we compute zd_lw(:,:2) and zd_up(:,:jpkm1). These values are
787      !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
788      !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
789      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
790         zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
791         zavt  = zsqen * zstt(ji,jj,jk)
792         zavm  = zsqen * zstm(ji,jj,jk)
793         p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
794         p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
795      END_3D
796      p_avt(A2D(nn_hls),1) = 0._wp
797      !
798      IF(sn_cfctl%l_prtctl) THEN
799         CALL prt_ctl( tab3d_1=en   , clinfo1=' gls  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
800         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls  - m: ', kdim=jpk )
801      ENDIF
802      !
803   END SUBROUTINE zdf_gls
804
805
806   SUBROUTINE zdf_gls_init
807      !!----------------------------------------------------------------------
808      !!                  ***  ROUTINE zdf_gls_init  ***
809      !!
810      !! ** Purpose :   Initialization of the vertical eddy diffivity and
811      !!              viscosity computed using a GLS turbulent closure scheme
812      !!
813      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
814      !!
815      !! ** input   :   Namlist namzdf_gls
816      !!
817      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
818      !!
819      !!----------------------------------------------------------------------
820      INTEGER ::   jk    ! dummy loop indices
821      INTEGER ::   ios   ! Local integer output status for namelist read
822      REAL(wp)::   zcr   ! local scalar
823      !!
824      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       &
825         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   &
826         &            nn_mxlice, rn_crban, rn_charn, rn_frac_hs,   &
827         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, &
828         &            nn_stab_func, nn_clos
829      !!----------------------------------------------------------
830      !
831      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
832901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_gls in reference namelist' )
833
834      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
835902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist' )
836      IF(lwm) WRITE ( numond, namzdf_gls )
837
838      IF(lwp) THEN                     !* Control print
839         WRITE(numout,*)
840         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
841         WRITE(numout,*) '~~~~~~~~~~~~'
842         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
843         WRITE(numout,*) '      minimum value of en                           rn_emin        = ', rn_emin
844         WRITE(numout,*) '      minimum value of eps                          rn_epsmin      = ', rn_epsmin
845         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim  = ', ln_length_lim
846         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp   = ', rn_clim_galp
847         WRITE(numout,*) '      TKE Surface boundary condition                nn_bc_surf     = ', nn_bc_surf
848         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
849         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
850         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
851         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
852         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
853         WRITE(numout,*) '      surface wave breaking under ice               nn_z0_ice      = ', nn_z0_ice
854         SELECT CASE( nn_z0_ice )
855         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on surface wave breaking'
856         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )'
857         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-fr_i(:,:)'
858         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )'
859         CASE DEFAULT
860            CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3')
861         END SELECT
862         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs
863         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
864         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
865         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
866         WRITE(numout,*) '      type of scaling under sea-ice                 nn_mxlice      = ', nn_mxlice
867         IF( nn_mxlice == 1 ) &
868            WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri        = ', rn_hsri
869         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
870            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice'
871            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness'
872            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean     sea-ice thickness'
873            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max      sea-ice thickness'
874            CASE DEFAULT
875               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 ')
876         END SELECT
877         WRITE(numout,*)
878      ENDIF
879
880      !                                !* allocate GLS arrays
881      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
882
883      !                                !* Check of some namelist values
884      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )              CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
885      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )              CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
886      IF( nn_z0_met  < 0   .OR. nn_z0_met    > 3 )              CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )
887      IF( nn_z0_met == 3  .AND. .NOT. (ln_wave .AND. ln_sdw ) ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T and ln_sdw=T' )
888      IF( nn_stab_func < 0 .OR. nn_stab_func > 3 )              CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )
889      IF( nn_clos      < 0 .OR. nn_clos      > 3 )              CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
890
891      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
892      !
893      CASE( 0 )                              ! k-kl  (Mellor-Yamada)
894         !
895         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
896         IF(lwp) WRITE(numout,*)
897         rpp     = 0._wp
898         rmm     = 1._wp
899         rnn     = 1._wp
900         rsc_tke = 1.96_wp
901         rsc_psi = 1.96_wp
902         rpsi1   = 0.9_wp
903         rpsi3p  = 1._wp
904         rpsi2   = 0.5_wp
905         !
906         SELECT CASE ( nn_stab_func )
907         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
908         CASE( 2 )      ;   rpsi3m = 2.62_wp       ! Canuto A stability functions
909         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
910         END SELECT
911         !
912      CASE( 1 )                              ! k-eps
913         !
914         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen'
915         IF(lwp) WRITE(numout,*)
916         rpp     =  3._wp
917         rmm     =  1.5_wp
918         rnn     = -1._wp
919         rsc_tke =  1._wp
920         rsc_psi =  1.2_wp  ! Schmidt number for psi
921         rpsi1   =  1.44_wp
922         rpsi3p  =  1._wp
923         rpsi2   =  1.92_wp
924         !
925         SELECT CASE ( nn_stab_func )
926         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
927         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
928         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
929         END SELECT
930         !
931      CASE( 2 )                              ! k-omega
932         !
933         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen'
934         IF(lwp) WRITE(numout,*)
935         rpp     = -1._wp
936         rmm     =  0.5_wp
937         rnn     = -1._wp
938         rsc_tke =  2._wp
939         rsc_psi =  2._wp
940         rpsi1   =  0.555_wp
941         rpsi3p  =  1._wp
942         rpsi2   =  0.833_wp
943         !
944         SELECT CASE ( nn_stab_func )
945         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
946         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
947         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
948         END SELECT
949         !
950      CASE( 3 )                              ! generic
951         !
952         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen'
953         IF(lwp) WRITE(numout,*)
954         rpp     = 2._wp
955         rmm     = 1._wp
956         rnn     = -0.67_wp
957         rsc_tke = 0.8_wp
958         rsc_psi = 1.07_wp
959         rpsi1   = 1._wp
960         rpsi3p  = 1._wp
961         rpsi2   = 1.22_wp
962         !
963         SELECT CASE ( nn_stab_func )
964         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
965         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
966         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
967         END SELECT
968         !
969      END SELECT
970
971      !
972      SELECT CASE ( nn_stab_func )     !* set the parameters of the stability functions
973      !
974      CASE ( 0 )                             ! Galperin stability functions
975         !
976         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin'
977         rc2     =  0._wp
978         rc3     =  0._wp
979         rc_diff =  1._wp
980         rc0     =  0.5544_wp
981         rcm_sf  =  0.9884_wp
982         rghmin  = -0.28_wp
983         rgh0    =  0.0233_wp
984         rghcri  =  0.02_wp
985         !
986      CASE ( 1 )                             ! Kantha-Clayson stability functions
987         !
988         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson'
989         rc2     =  0.7_wp
990         rc3     =  0.2_wp
991         rc_diff =  1._wp
992         rc0     =  0.5544_wp
993         rcm_sf  =  0.9884_wp
994         rghmin  = -0.28_wp
995         rgh0    =  0.0233_wp
996         rghcri  =  0.02_wp
997         !
998      CASE ( 2 )                             ! Canuto A stability functions
999         !
1000         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A'
1001         rs0 = 1.5_wp * rl1 * rl5*rl5
1002         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
1003         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
1004         rs4 = 2._wp * rl5
1005         rs5 = 2._wp * rl4
1006         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
1007            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
1008         rd0 = 3._wp * rl5*rl5
1009         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
1010         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
1011         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
1012         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
1013         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
1014         rc0 = 0.5268_wp
1015         rf6 = 8._wp / (rc0**6._wp)
1016         rc_diff = SQRT(2._wp) / (rc0**3._wp)
1017         rcm_sf  =  0.7310_wp
1018         rghmin  = -0.28_wp
1019         rgh0    =  0.0329_wp
1020         rghcri  =  0.03_wp
1021         !
1022      CASE ( 3 )                             ! Canuto B stability functions
1023         !
1024         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B'
1025         rs0 = 1.5_wp * rm1 * rm5*rm5
1026         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
1027         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
1028         rs4 = 2._wp * rm5
1029         rs5 = 2._wp * rm4
1030         rs6 = (2._wp/3._wp) * rm5 * (3._wp*rm3*rm3-rm2*rm2) - 0.5_wp * rm5*rm1*(3._wp*rm3-rm2) + 0.75_wp * rm1*(rm6-rm7)
1031         rd0 = 3._wp * rm5*rm5
1032         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
1033         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
1034         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
1035         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
1036         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
1037         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
1038         rf6 = 8._wp / ( rc0**6._wp )
1039         rc_diff = SQRT(2._wp)/(rc0**3.)
1040         rcm_sf  =  0.7470_wp
1041         rghmin  = -0.28_wp
1042         rgh0    =  0.0444_wp
1043         rghcri  =  0.0414_wp
1044         !
1045      END SELECT
1046
1047      !                                !* Set Schmidt number for psi diffusion in the wave breaking case
1048      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1049      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1050      IF( ln_sigpsi ) THEN
1051         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
1052         ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
1053         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
1054         ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn
1055         rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.))
1056      ELSE
1057         rsc_psi0 = rsc_psi
1058      ENDIF
1059
1060      !                                !* Shear free turbulence parameters
1061      !
1062      ra_sf  = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) &
1063               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
1064
1065      IF ( rn_crban==0._wp ) THEN
1066         rl_sf = vkarmn
1067      ELSE
1068         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp) * rsc_tke        &
1069            &                                            + 12._wp*rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm)   &
1070            &                                                     *SQRT(rsc_tke*(rsc_tke                 &
1071            &                                                        + 24._wp*rsc_psi0*rpsi2)) )         &
1072            &                                              /(12._wp*rnn**2.)                             )
1073      ENDIF
1074
1075      !
1076      IF(lwp) THEN                     !* Control print
1077         WRITE(numout,*)
1078         WRITE(numout,*) '   Limit values :'
1079         WRITE(numout,*) '      Parameter  m = ', rmm
1080         WRITE(numout,*) '      Parameter  n = ', rnn
1081         WRITE(numout,*) '      Parameter  p = ', rpp
1082         WRITE(numout,*) '      rpsi1    = ', rpsi1
1083         WRITE(numout,*) '      rpsi2    = ', rpsi2
1084         WRITE(numout,*) '      rpsi3m   = ', rpsi3m
1085         WRITE(numout,*) '      rpsi3p   = ', rpsi3p
1086         WRITE(numout,*) '      rsc_tke  = ', rsc_tke
1087         WRITE(numout,*) '      rsc_psi  = ', rsc_psi
1088         WRITE(numout,*) '      rsc_psi0 = ', rsc_psi0
1089         WRITE(numout,*) '      rc0      = ', rc0
1090         WRITE(numout,*)
1091         WRITE(numout,*) '   Shear free turbulence parameters:'
1092         WRITE(numout,*) '      rcm_sf   = ', rcm_sf
1093         WRITE(numout,*) '      ra_sf    = ', ra_sf
1094         WRITE(numout,*) '      rl_sf    = ', rl_sf
1095      ENDIF
1096
1097      !                                !* Constants initialization
1098      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
1099      rc03  = rc02 * rc0
1100      rc04  = rc03 * rc0
1101      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
1102      rsbc_tke2 = rn_Dt * rn_crban / rl_sf                                 ! Neumann + Wave breaking
1103      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
1104      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer
1105      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
1106      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
1107      rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi
1108      rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
1109      !
1110      rfact_tke = -0.5_wp / rsc_tke * rn_Dt                                ! Cst used for the Diffusion term of tke
1111      rfact_psi = -0.5_wp / rsc_psi * rn_Dt                                ! Cst used for the Diffusion term of tke
1112      !
1113      !                                !* Wall proximity function
1114!!gm tmask or wmask ????
1115      zwall(:,:,:) = 1._wp * tmask(:,:,:)
1116
1117      !                                !* read or initialize all required files
1118      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxl_n)
1119      !
1120   END SUBROUTINE zdf_gls_init
1121
1122
1123   SUBROUTINE gls_rst( kt, cdrw )
1124      !!---------------------------------------------------------------------
1125      !!                   ***  ROUTINE gls_rst  ***
1126      !!
1127      !! ** Purpose :   Read or write TKE file (en) in restart file
1128      !!
1129      !! ** Method  :   use of IOM library
1130      !!                if the restart does not contain TKE, en is either
1131      !!                set to rn_emin or recomputed (nn_igls/=0)
1132      !!----------------------------------------------------------------------
1133      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
1134      !!
1135      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1136      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1137      !
1138      INTEGER ::   jit, jk   ! dummy loop indices
1139      INTEGER ::   id1, id2, id3, id4
1140      INTEGER ::   ji, jj, ikbu, ikbv
1141      REAL(wp)::   cbx, cby
1142      !!----------------------------------------------------------------------
1143      !
1144      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1145         !                                   ! ---------------
1146         IF( ln_rstart ) THEN                   !* Read the restart file
1147            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. )
1148            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
1149            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
1150            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
1151            !
1152            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
1153               CALL iom_get( numror, jpdom_auto, 'en'    , en     )
1154               CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k  )
1155               CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k  )
1156               CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n )
1157            ELSE
1158               IF(lwp) WRITE(numout,*)
1159               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values'
1160               en    (:,:,:) = rn_emin
1161               hmxl_n(:,:,:) = 0.05_wp
1162               ! avt_k, avm_k already set to the background value in zdf_phy_init
1163            ENDIF
1164         ELSE                                   !* Start from rest
1165            IF(lwp) WRITE(numout,*)
1166            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values'
1167            en    (:,:,:) = rn_emin
1168            hmxl_n(:,:,:) = 0.05_wp
1169            ! avt_k, avm_k already set to the background value in zdf_phy_init
1170         ENDIF
1171         !
1172      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1173         !                                   ! -------------------
1174         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1175         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en     )
1176         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k  )
1177         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k  )
1178         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
1179         !
1180      ENDIF
1181      !
1182   END SUBROUTINE gls_rst
1183
1184   !!======================================================================
1185END MODULE zdfgls
Note: See TracBrowser for help on using the repository browser.