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 @ 15071

Last change on this file since 15071 was 15071, checked in by clem, 3 years ago

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