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.
zdfgls.F90 in NEMO/trunk/src/OCE/ZDF – NEMO

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

Last change on this file since 14983 was 14967, checked in by cbricaud, 3 years ago

introduce limit conditions under ice caps in GLS; see #2593 and #2604
  • Property svn:keywords set to Id
File size: 65.8 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         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
357            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
358               IF( mikt(ji,jj) > 1 )THEN
359                  itop   = mikt(ji,jj)       ! k   top w-point
360                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
361                  !                                                ! mask at the
362                  !                                                ocean surface
363                  !                                                points
364                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
365                  !
366                  ! Dirichlet condition applied at:
367                  !     top level (itop)         &      Just below it (itopp1)
368                  zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
369                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
370                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
371                  en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
372               ENDIF
373            END_2D
374         ENDIF
375         !
376      CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
377         !
378         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
379            ! Dirichlet conditions at k=1
380            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  )
381            zd_lw(ji,jj,1) = en(ji,jj,1)
382            zd_up(ji,jj,1) = 0._wp
383            zdiag(ji,jj,1) = 1._wp
384            !
385            ! at k=2, set de/dz=Fw
386            !cbr
387            ! zdiag zd_lw not defined/used on the halo
388            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
389            zd_lw(ji,jj,2) = 0._wp
390            !
391            zkar (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) ))
392            zflxs(ji,jj)   = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) &
393                &                    * (  ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj)  )**(1.5_wp*ra_sf)
394!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf)
395            en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
396         END_2D
397         !
398         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
399            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
400               IF( mikt(ji,jj) > 1 )THEN
401                  itop   = mikt(ji,jj)       ! k   top w-point
402                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
403                  !                                                ! mask at the
404                  !                                                ocean surface
405                  !                                                points
406                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
407                  !
408                  ! Bottom level Dirichlet condition:
409                  !     Bottom level (ibot)      &      Just above it (ibotm1)
410                  !         Dirichlet            !         Neumann
411                  zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
412                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
413                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
414                  en   (ji,jj,itop) = z_en
415               ENDIF
416            END_2D
417         ENDIF
418         !
419      END SELECT
420
421      ! Bottom boundary condition on tke
422      ! --------------------------------
423      !
424      SELECT CASE ( nn_bc_bot )
425      !
426      CASE ( 0 )             ! Dirichlet
427         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
428         !                      ! Balance between the production and the dissipation terms
429         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
430!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
431!!   With thick deep ocean level thickness, this may be quite large, no ???
432!!   in particular in ocean cavities where top stratification can be large...
433            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
434            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
435            !
436            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
437            !
438            ! Dirichlet condition applied at:
439            !     Bottom level (ibot)      &      Just above it (ibotm1)
440            zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp
441            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
442            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp
443            en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en
444         END_2D
445         !
446         ! NOTE: ctl_stop with ln_isfcav when using GLS
447         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
448            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
449               itop   = mikt(ji,jj)       ! k   top w-point
450               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
451               !                                                ! mask at the ocean surface points
452               z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
453               !
454 !!gm TO BE VERIFIED !!!
455               ! Dirichlet condition applied at:
456               !     top level (itop)         &      Just below it (itopp1)
457               zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
458               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
459               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
460               en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
461            END_2D
462         ENDIF
463         !
464      CASE ( 1 )             ! Neumman boundary condition
465         !
466         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
467            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
468            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
469            !
470            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
471            !
472            ! Bottom level Dirichlet condition:
473            !     Bottom level (ibot)      &      Just above it (ibotm1)
474            !         Dirichlet            !         Neumann
475            zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag
476            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
477            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
478            en   (ji,jj,ibot) = z_en
479         END_2D
480         ! NOTE: ctl_stop with ln_isfcav when using GLS
481         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
482            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
483               itop   = mikt(ji,jj)       ! k   top w-point
484               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
485               !                                                ! mask at the ocean surface points
486               z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
487               !
488               ! Bottom level Dirichlet condition:
489               !     Bottom level (ibot)      &      Just above it (ibotm1)
490               !         Dirichlet            !         Neumann
491               zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
492               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
493               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
494               en   (ji,jj,itop) = z_en
495            END_2D
496         ENDIF
497         !
498      END SELECT
499
500      ! Matrix inversion (en prescribed at surface and the bottom)
501      ! ----------------------------------------------------------
502      !
503      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
504         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
505      END_3D
506      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
507         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)
508      END_3D
509      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
510         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
511      END_3D
512      !                                            ! set the minimum value of tke
513      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
514         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin )
515      END_3D
516
517      !!----------------------------------------!!
518      !!   Solve prognostic equation for psi    !!
519      !!----------------------------------------!!
520
521      ! Set psi to previous time step value
522      !
523      SELECT CASE ( nn_clos )
524      !
525      CASE( 0 )               ! k-kl  (Mellor-Yamada)
526         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
527            psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
528         END_3D
529         !
530      CASE( 1 )               ! k-eps
531         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
532            psi(ji,jj,jk)  = eps(ji,jj,jk)
533         END_3D
534         !
535      CASE( 2 )               ! k-w
536         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
537            psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
538         END_3D
539         !
540      CASE( 3 )               ! generic
541         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
542            psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn
543         END_3D
544         !
545      END SELECT
546      !
547      ! Now gls (output in psi)
548      ! -------------------------------
549      ! Resolution of a tridiagonal linear system by a "methode de chasse"
550      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
551      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
552      ! Warning : after this step, en : right hand side of the matrix
553
554      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
555         !
556         ! psi / k
557         zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
558         !
559         ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable)
560         zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
561         !
562         rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p
563         !
564         ! shear prod. - stratif. destruction
565         prod = rpsi1 * zratio * p_sh2(ji,jj,jk)
566         !
567         ! stratif. destruction
568         buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
569         !
570         ! shear prod. - stratif. destruction
571         diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
572         !
573         zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
574         !
575         zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term
576         zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
577         !
578         ! building the matrix
579         zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
580         !                                               ! lower diagonal
581         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) )   &
582            &            / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
583         !                                               ! upper diagonal
584         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) )   &
585            &            / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
586         !                                               ! diagonal
587         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk)
588         !                                               ! right hand side in psi
589         psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk)
590      END_3D
591      !
592      zdiag(:,:,jpk) = 1._wp
593
594      ! Surface boundary condition on psi
595      ! ---------------------------------
596      !
597      SELECT CASE ( nn_bc_surf )
598      !
599      CASE ( 0 )             ! Dirichlet boundary conditions
600         !
601         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
602            ! Surface value
603            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf ! Cosmetic
604            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
605            zd_lw(ji,jj,1) = psi(ji,jj,1)
606            zd_up(ji,jj,1) = 0._wp
607            zdiag(ji,jj,1) = 1._wp
608            !
609            ! One level below
610            zkar    (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) )))
611            zdep    (ji,jj)   = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj)
612            psi     (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
613            zd_lw(ji,jj,2) = 0._wp
614            zd_up(ji,jj,2) = 0._wp
615            zdiag(ji,jj,2) = 1._wp
616         END_2D
617         !
618      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
619         !
620         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
621            ! Surface value: Dirichlet
622            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf
623            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
624            zd_lw(ji,jj,1) = psi(ji,jj,1)
625            zd_up(ji,jj,1) = 0._wp
626            zdiag(ji,jj,1) = 1._wp
627            !
628            ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo
629            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
630            zd_lw(ji,jj,2) = 0._wp
631            !
632            ! Set psi vertical flux at the surface:
633            zkar (ji,jj)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope
634            zdep (ji,jj)   = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf)
635            zflxs(ji,jj)   = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) &
636               &           *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp)
637            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)) * &
638               &           ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.)
639            zflxs(ji,jj)   = zdep(ji,jj) * zflxs(ji,jj)
640            psi  (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
641         END_2D
642         !
643      END SELECT
644
645      ! Bottom boundary condition on psi
646      ! --------------------------------
647      !
648!!gm should be done for ISF (top boundary cond.)
649!!gm so, totally new staff needed      ===>>> think about that !
650!
651      SELECT CASE ( nn_bc_bot )     ! bottom boundary
652      !
653      CASE ( 0 )             ! Dirichlet
654         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
655         !                      ! Balance between the production and the dissipation terms
656         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
657            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
658            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
659            zdep(ji,jj) = vkarmn * r_z0_bot
660            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
661            zd_lw(ji,jj,ibot) = 0._wp
662            zd_up(ji,jj,ibot) = 0._wp
663            zdiag(ji,jj,ibot) = 1._wp
664            !
665            ! Just above last level, Dirichlet condition again (GOTM like)
666            zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) )
667            psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
668            zd_lw(ji,jj,ibotm1) = 0._wp
669            zd_up(ji,jj,ibotm1) = 0._wp
670            zdiag(ji,jj,ibotm1) = 1._wp
671         END_2D
672         !
673         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
674            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
675               IF ( mikt(ji,jj) > 1 ) THEN
676                  itop   = mikt(ji,jj)       ! k   top w-point
677                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
678                  !
679                  zdep(ji,jj) = vkarmn * r_z0_top
680                  psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn
681                  zd_lw(ji,jj,itop) = 0._wp
682                  zd_up(ji,jj,itop) = 0._wp
683                  zdiag(ji,jj,itop) = 1._wp
684                  !
685                  ! Just above last level, Dirichlet condition again (GOTM like)
686                  zdep(ji,jj) = vkarmn * ( r_z0_top + e3t(ji,jj,itopp1,Kmm) )
687                  psi (ji,jj,itopp1) = rc0**rpp * en(ji,jj,itop  )**rmm *zdep(ji,jj)**rnn
688                  zd_lw(ji,jj,itopp1) = 0._wp
689                  zd_up(ji,jj,itopp1) = 0._wp
690                  zdiag(ji,jj,itopp1) = 1._wp
691               END IF
692            END_2D
693         END IF
694         !
695      CASE ( 1 )             ! Neumman boundary condition
696         !
697         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
698            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
699            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
700            !
701            ! Bottom level Dirichlet condition:
702            zdep(ji,jj) = vkarmn * r_z0_bot
703            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
704            !
705            zd_lw(ji,jj,ibot) = 0._wp
706            zd_up(ji,jj,ibot) = 0._wp
707            zdiag(ji,jj,ibot) = 1._wp
708            !
709            ! Just above last level: Neumann condition with flux injection
710            zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag
711            zd_up(ji,jj,ibotm1) = 0.
712            !
713            ! Set psi vertical flux at the bottom:
714            zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm)
715            zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
716               &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
717            psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm)
718         END_2D
719         !
720         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
721            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
722               IF ( mikt(ji,jj) > 1 ) THEN
723                  itop   = mikt(ji,jj)       ! k   top w-point
724                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
725                  !
726                  ! Bottom level Dirichlet condition:
727                  zdep(ji,jj) = vkarmn * r_z0_top
728                  psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn
729                  !
730                  zd_lw(ji,jj,itop) = 0._wp
731                  zd_up(ji,jj,itop) = 0._wp
732                  zdiag(ji,jj,itop) = 1._wp
733                  !
734                  ! Just below cavity level: Neumann condition with flux
735                  ! injection
736                  zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) ! Remove zd_up from zdiag
737                  zd_up(ji,jj,itopp1) = 0._wp
738                  !
739                  ! Set psi vertical flux below cavity:
740                  zdep(ji,jj) = r_z0_top + 0.5_wp*e3t(ji,jj,itopp1,Kmm)
741                  zflxb = rsbc_psi2 * ( p_avm(ji,jj,itop) + p_avm(ji,jj,itopp1))   &
742                     &  * (0.5_wp*(en(ji,jj,itop)+en(ji,jj,itopp1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
743                  psi(ji,jj,itopp1) = psi(ji,jj,itopp1) + zflxb / e3w(ji,jj,itopp1,Kmm)
744               END IF
745            END_2D
746         END IF
747
748         !
749      END SELECT
750
751      ! Matrix inversion
752      ! ----------------
753      !
754      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
755         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
756      END_3D
757      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
758         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)
759      END_3D
760      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
761         psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
762      END_3D
763
764      ! Set dissipation
765      !----------------
766
767      SELECT CASE ( nn_clos )
768      !
769      CASE( 0 )               ! k-kl  (Mellor-Yamada)
770         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
771            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)
772         END_3D
773         !
774      CASE( 1 )               ! k-eps
775         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
776            eps(ji,jj,jk) = psi(ji,jj,jk)
777         END_3D
778         !
779      CASE( 2 )               ! k-w
780         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
781            eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
782         END_3D
783         !
784      CASE( 3 )               ! generic
785         zcoef = rc0**( 3._wp  + rpp/rnn )
786         zex1  =      ( 1.5_wp + rmm/rnn )
787         zex2  = -1._wp / rnn
788         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
789            eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
790         END_3D
791         !
792      END SELECT
793
794      ! Limit dissipation rate under stable stratification
795      ! --------------------------------------------------
796      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
797         ! limitation
798         eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
799         hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
800      END_3D
801      IF( ln_length_lim ) THEN        ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
802         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
803            zrn2 = MAX( rn2(ji,jj,jk), rsmall )
804            hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
805         END_3D
806      ENDIF
807
808      !
809      ! Stability function and vertical viscosity and diffusivity
810      ! ---------------------------------------------------------
811      !
812      SELECT CASE ( nn_stab_func )
813      !
814      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
815         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
816            ! zcof =  l²/q²
817            zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
818            ! Gh = -N²l²/q²
819            gh = - rn2(ji,jj,jk) * zcof
820            gh = MIN( gh, rgh0   )
821            gh = MAX( gh, rghmin )
822            ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
823            sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) )
824            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)
825            !
826            ! Store stability function in zstt and zstm
827            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
828            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
829         END_3D
830         !
831      CASE ( 2, 3 )               ! Canuto stability functions
832         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
833            ! zcof =  l²/q²
834            zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
835            ! Gh = -N²l²/q²
836            gh = - rn2(ji,jj,jk) * zcof
837            gh = MIN( gh, rgh0   )
838            gh = MAX( gh, rghmin )
839            gh = gh * rf6
840            ! Gm =  M²l²/q² Shear number
841            shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
842            gm = MAX( shr * zcof , 1.e-10 )
843            gm = gm * rf6
844            gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm )
845            ! Stability functions from Canuto
846            rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm
847            sm = (rs0 - rs1*gh + rs2*gm) / rcff
848            sh = (rs4 - rs5*gh + rs6*gm) / rcff
849            !
850            ! Store stability function in zstt and zstm
851            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
852            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
853         END_3D
854         !
855      END SELECT
856
857      ! Boundary conditions on stability functions for momentum (Neumann):
858      ! Lines below are useless if GOTM style Dirichlet conditions are used
859
860      zstm(:,:,1) = zstm(:,:,2)
861
862      ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
863      zstm(:,:,jpk) = 0.
864      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )             ! update bottom with good values
865         zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
866      END_2D
867
868      zstt(:,:,  1) = wmask(A2D(nn_hls),  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
869      zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
870
871!!gm should be done for ISF (top boundary cond.)
872!!gm so, totally new staff needed!!gm
873
874      ! Compute diffusivities/viscosities
875      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
876      !  -> yes BUT p_avm(:,:1) and p_avm(:,:jpk) are used when we compute zd_lw(:,:2) and zd_up(:,:jpkm1). These values are
877      !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
878      !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
879      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
880         zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
881         zavt  = zsqen * zstt(ji,jj,jk)
882         zavm  = zsqen * zstm(ji,jj,jk)
883         p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
884         p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
885      END_3D
886      p_avt(A2D(nn_hls),1) = 0._wp
887      !
888      IF(sn_cfctl%l_prtctl) THEN
889         CALL prt_ctl( tab3d_1=en   , clinfo1=' gls  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
890         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls  - m: ', kdim=jpk )
891      ENDIF
892      !
893   END SUBROUTINE zdf_gls
894
895
896   SUBROUTINE zdf_gls_init
897      !!----------------------------------------------------------------------
898      !!                  ***  ROUTINE zdf_gls_init  ***
899      !!
900      !! ** Purpose :   Initialization of the vertical eddy diffivity and
901      !!              viscosity computed using a GLS turbulent closure scheme
902      !!
903      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
904      !!
905      !! ** input   :   Namlist namzdf_gls
906      !!
907      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
908      !!
909      !!----------------------------------------------------------------------
910      INTEGER ::   jk    ! dummy loop indices
911      INTEGER ::   ios   ! Local integer output status for namelist read
912      REAL(wp)::   zcr   ! local scalar
913      !!
914      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       &
915         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   &
916         &            nn_mxlice, rn_crban, rn_charn, rn_frac_hs,   &
917         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, &
918         &            nn_stab_func, nn_clos
919      !!----------------------------------------------------------
920      !
921      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
922901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_gls in reference namelist' )
923
924      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
925902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist' )
926      IF(lwm) WRITE ( numond, namzdf_gls )
927
928      IF(lwp) THEN                     !* Control print
929         WRITE(numout,*)
930         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
931         WRITE(numout,*) '~~~~~~~~~~~~'
932         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
933         WRITE(numout,*) '      minimum value of en                           rn_emin        = ', rn_emin
934         WRITE(numout,*) '      minimum value of eps                          rn_epsmin      = ', rn_epsmin
935         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim  = ', ln_length_lim
936         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp   = ', rn_clim_galp
937         WRITE(numout,*) '      TKE Surface boundary condition                nn_bc_surf     = ', nn_bc_surf
938         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
939         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
940         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
941         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
942         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
943         WRITE(numout,*) '      surface wave breaking under ice               nn_z0_ice      = ', nn_z0_ice
944         SELECT CASE( nn_z0_ice )
945         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on surface wave breaking'
946         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )'
947         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-fr_i(:,:)'
948         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )'
949         CASE DEFAULT
950            CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3')
951         END SELECT
952         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs
953         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
954         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
955         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
956         WRITE(numout,*) '      type of scaling under sea-ice                 nn_mxlice      = ', nn_mxlice
957         IF( nn_mxlice == 1 ) &
958            WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri        = ', rn_hsri
959         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
960            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice'
961            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness'
962            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean     sea-ice thickness'
963            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max      sea-ice thickness'
964            CASE DEFAULT
965               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 ')
966         END SELECT
967         WRITE(numout,*)
968      ENDIF
969
970      !                                !* allocate GLS arrays
971      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
972
973      !                                !* Check of some namelist values
974      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' )
975      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' )
976      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' )
977      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' )
978      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' )
979      IF( nn_clos      < 0 .OR. nn_clos      > 3 )              CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
980
981      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
982      !
983      CASE( 0 )                              ! k-kl  (Mellor-Yamada)
984         !
985         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
986         IF(lwp) WRITE(numout,*)
987         rpp     = 0._wp
988         rmm     = 1._wp
989         rnn     = 1._wp
990         rsc_tke = 1.96_wp
991         rsc_psi = 1.96_wp
992         rpsi1   = 0.9_wp
993         rpsi3p  = 1._wp
994         rpsi2   = 0.5_wp
995         !
996         SELECT CASE ( nn_stab_func )
997         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
998         CASE( 2 )      ;   rpsi3m = 2.62_wp       ! Canuto A stability functions
999         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
1000         END SELECT
1001         !
1002      CASE( 1 )                              ! k-eps
1003         !
1004         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen'
1005         IF(lwp) WRITE(numout,*)
1006         rpp     =  3._wp
1007         rmm     =  1.5_wp
1008         rnn     = -1._wp
1009         rsc_tke =  1._wp
1010         rsc_psi =  1.2_wp  ! Schmidt number for psi
1011         rpsi1   =  1.44_wp
1012         rpsi3p  =  1._wp
1013         rpsi2   =  1.92_wp
1014         !
1015         SELECT CASE ( nn_stab_func )
1016         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
1017         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
1018         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
1019         END SELECT
1020         !
1021      CASE( 2 )                              ! k-omega
1022         !
1023         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen'
1024         IF(lwp) WRITE(numout,*)
1025         rpp     = -1._wp
1026         rmm     =  0.5_wp
1027         rnn     = -1._wp
1028         rsc_tke =  2._wp
1029         rsc_psi =  2._wp
1030         rpsi1   =  0.555_wp
1031         rpsi3p  =  1._wp
1032         rpsi2   =  0.833_wp
1033         !
1034         SELECT CASE ( nn_stab_func )
1035         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
1036         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
1037         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
1038         END SELECT
1039         !
1040      CASE( 3 )                              ! generic
1041         !
1042         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen'
1043         IF(lwp) WRITE(numout,*)
1044         rpp     = 2._wp
1045         rmm     = 1._wp
1046         rnn     = -0.67_wp
1047         rsc_tke = 0.8_wp
1048         rsc_psi = 1.07_wp
1049         rpsi1   = 1._wp
1050         rpsi3p  = 1._wp
1051         rpsi2   = 1.22_wp
1052         !
1053         SELECT CASE ( nn_stab_func )
1054         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
1055         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
1056         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
1057         END SELECT
1058         !
1059      END SELECT
1060
1061      !
1062      SELECT CASE ( nn_stab_func )     !* set the parameters of the stability functions
1063      !
1064      CASE ( 0 )                             ! Galperin stability functions
1065         !
1066         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin'
1067         rc2     =  0._wp
1068         rc3     =  0._wp
1069         rc_diff =  1._wp
1070         rc0     =  0.5544_wp
1071         rcm_sf  =  0.9884_wp
1072         rghmin  = -0.28_wp
1073         rgh0    =  0.0233_wp
1074         rghcri  =  0.02_wp
1075         !
1076      CASE ( 1 )                             ! Kantha-Clayson stability functions
1077         !
1078         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson'
1079         rc2     =  0.7_wp
1080         rc3     =  0.2_wp
1081         rc_diff =  1._wp
1082         rc0     =  0.5544_wp
1083         rcm_sf  =  0.9884_wp
1084         rghmin  = -0.28_wp
1085         rgh0    =  0.0233_wp
1086         rghcri  =  0.02_wp
1087         !
1088      CASE ( 2 )                             ! Canuto A stability functions
1089         !
1090         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A'
1091         rs0 = 1.5_wp * rl1 * rl5*rl5
1092         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
1093         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
1094         rs4 = 2._wp * rl5
1095         rs5 = 2._wp * rl4
1096         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
1097            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
1098         rd0 = 3._wp * rl5*rl5
1099         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
1100         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
1101         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
1102         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
1103         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
1104         rc0 = 0.5268_wp
1105         rf6 = 8._wp / (rc0**6._wp)
1106         rc_diff = SQRT(2._wp) / (rc0**3._wp)
1107         rcm_sf  =  0.7310_wp
1108         rghmin  = -0.28_wp
1109         rgh0    =  0.0329_wp
1110         rghcri  =  0.03_wp
1111         !
1112      CASE ( 3 )                             ! Canuto B stability functions
1113         !
1114         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B'
1115         rs0 = 1.5_wp * rm1 * rm5*rm5
1116         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
1117         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
1118         rs4 = 2._wp * rm5
1119         rs5 = 2._wp * rm4
1120         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)
1121         rd0 = 3._wp * rm5*rm5
1122         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
1123         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
1124         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
1125         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
1126         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
1127         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
1128         rf6 = 8._wp / ( rc0**6._wp )
1129         rc_diff = SQRT(2._wp)/(rc0**3.)
1130         rcm_sf  =  0.7470_wp
1131         rghmin  = -0.28_wp
1132         rgh0    =  0.0444_wp
1133         rghcri  =  0.0414_wp
1134         !
1135      END SELECT
1136
1137      !                                !* Set Schmidt number for psi diffusion in the wave breaking case
1138      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1139      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1140      IF( ln_sigpsi ) THEN
1141         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
1142         ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
1143         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
1144         ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn
1145         rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.))
1146      ELSE
1147         rsc_psi0 = rsc_psi
1148      ENDIF
1149
1150      !                                !* Shear free turbulence parameters
1151      !
1152      ra_sf  = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) &
1153               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
1154
1155      IF ( rn_crban==0._wp ) THEN
1156         rl_sf = vkarmn
1157      ELSE
1158         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp) * rsc_tke        &
1159            &                                            + 12._wp*rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm)   &
1160            &                                                     *SQRT(rsc_tke*(rsc_tke                 &
1161            &                                                        + 24._wp*rsc_psi0*rpsi2)) )         &
1162            &                                              /(12._wp*rnn**2.)                             )
1163      ENDIF
1164
1165      !
1166      IF(lwp) THEN                     !* Control print
1167         WRITE(numout,*)
1168         WRITE(numout,*) '   Limit values :'
1169         WRITE(numout,*) '      Parameter  m = ', rmm
1170         WRITE(numout,*) '      Parameter  n = ', rnn
1171         WRITE(numout,*) '      Parameter  p = ', rpp
1172         WRITE(numout,*) '      rpsi1    = ', rpsi1
1173         WRITE(numout,*) '      rpsi2    = ', rpsi2
1174         WRITE(numout,*) '      rpsi3m   = ', rpsi3m
1175         WRITE(numout,*) '      rpsi3p   = ', rpsi3p
1176         WRITE(numout,*) '      rsc_tke  = ', rsc_tke
1177         WRITE(numout,*) '      rsc_psi  = ', rsc_psi
1178         WRITE(numout,*) '      rsc_psi0 = ', rsc_psi0
1179         WRITE(numout,*) '      rc0      = ', rc0
1180         WRITE(numout,*)
1181         WRITE(numout,*) '   Shear free turbulence parameters:'
1182         WRITE(numout,*) '      rcm_sf   = ', rcm_sf
1183         WRITE(numout,*) '      ra_sf    = ', ra_sf
1184         WRITE(numout,*) '      rl_sf    = ', rl_sf
1185      ENDIF
1186
1187      !                                !* Constants initialization
1188      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
1189      rc03  = rc02 * rc0
1190      rc04  = rc03 * rc0
1191      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
1192      rsbc_tke2 = rn_Dt * rn_crban / rl_sf                                 ! Neumann + Wave breaking
1193      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
1194      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer
1195      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
1196      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
1197      rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi
1198      rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
1199      !
1200      rfact_tke = -0.5_wp / rsc_tke * rn_Dt                                ! Cst used for the Diffusion term of tke
1201      rfact_psi = -0.5_wp / rsc_psi * rn_Dt                                ! Cst used for the Diffusion term of tke
1202      !
1203      !                                !* Wall proximity function
1204!!gm tmask or wmask ????
1205      zwall(:,:,:) = 1._wp * tmask(:,:,:)
1206
1207      !                                !* read or initialize all required files
1208      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxl_n)
1209      !
1210   END SUBROUTINE zdf_gls_init
1211
1212
1213   SUBROUTINE gls_rst( kt, cdrw )
1214      !!---------------------------------------------------------------------
1215      !!                   ***  ROUTINE gls_rst  ***
1216      !!
1217      !! ** Purpose :   Read or write TKE file (en) in restart file
1218      !!
1219      !! ** Method  :   use of IOM library
1220      !!                if the restart does not contain TKE, en is either
1221      !!                set to rn_emin or recomputed (nn_igls/=0)
1222      !!----------------------------------------------------------------------
1223      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
1224      !!
1225      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1226      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1227      !
1228      INTEGER ::   jit, jk   ! dummy loop indices
1229      INTEGER ::   id1, id2, id3, id4
1230      INTEGER ::   ji, jj, ikbu, ikbv
1231      REAL(wp)::   cbx, cby
1232      !!----------------------------------------------------------------------
1233      !
1234      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1235         !                                   ! ---------------
1236         IF( ln_rstart ) THEN                   !* Read the restart file
1237            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. )
1238            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
1239            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
1240            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
1241            !
1242            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
1243               CALL iom_get( numror, jpdom_auto, 'en'    , en     )
1244               CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k  )
1245               CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k  )
1246               CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n )
1247            ELSE
1248               IF(lwp) WRITE(numout,*)
1249               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values'
1250               en    (:,:,:) = rn_emin
1251               hmxl_n(:,:,:) = 0.05_wp
1252               ! avt_k, avm_k already set to the background value in zdf_phy_init
1253            ENDIF
1254         ELSE                                   !* Start from rest
1255            IF(lwp) WRITE(numout,*)
1256            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values'
1257            en    (:,:,:) = rn_emin
1258            hmxl_n(:,:,:) = 0.05_wp
1259            ! avt_k, avm_k already set to the background value in zdf_phy_init
1260         ENDIF
1261         !
1262      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1263         !                                   ! -------------------
1264         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1265         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en     )
1266         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k  )
1267         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k  )
1268         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
1269         !
1270      ENDIF
1271      !
1272   END SUBROUTINE gls_rst
1273
1274   !!======================================================================
1275END MODULE zdfgls
Note: See TracBrowser for help on using the repository browser.