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

Last change on this file since 14903 was 14834, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling

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