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 branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 7905

Last change on this file since 7905 was 7905, checked in by jcastill, 8 years ago

Series of small bug fixes and stetic changes:

-Fix possible bug in the calculation of Stokes-Coriolis
-Move all the wave control variables to namelist namsbc_wave
-Use one namelist variable instead of two to set Stokes drift velocity coupling
-Cap the values of the Craig and Banner constant as calculated from wave input fields to take into account small values of the friction velocity
-Add new Phillips parametrization for Stokes drift vertical velocity, using the inverse depth scale as in Breivik 2015, instead of the peak wave number as calculated from wave input fields
-Better control of the wave fields that are read from file depending on the wave parameters

File size: 67.7 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   !!----------------------------------------------------------------------
10#if defined key_zdfgls   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_zdfgls'                 Generic Length Scale vertical physics
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 zdf_oce        ! ocean vertical physics
22   USE zdfbfr         ! bottom friction (only for rn_bfrz0)
23   USE sbc_oce        ! surface boundary condition: ocean
24   USE phycst         ! physical constants
25   USE zdfmxl         ! mixed layer
26   USE sbcwave
27   !
28   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
29   USE lib_mpp        ! MPP manager
30   USE wrk_nemo       ! work arrays
31   USE prtctl         ! Print control
32   USE in_out_manager ! I/O manager
33   USE iom            ! I/O manager library
34   USE timing         ! Timing
35   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   zdf_gls        ! routine called in step module
41   PUBLIC   zdf_gls_init   ! routine called in opa module
42   PUBLIC   gls_rst        ! routine called in step module
43
44   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfgls = .TRUE.   !: TKE vertical mixing flag
45   !
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points
50
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_tke1
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_tke3
53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_psi1
54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rtrans
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_bc_surf        ! surface boundary condition (=0/1)
60   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1)
61   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation
62   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
63   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
64   REAL(wp) ::   rn_clim_galp      ! Holt 2008 value for k-eps: 0.267
65   REAL(wp) ::   rn_epsmin         ! minimum value of dissipation (m2/s3)
66   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
67   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
68   REAL(wp) ::   rn_crban_default  ! Craig and Banner constant for surface breaking waves mixing
69   REAL(wp) ::   rn_hsro           ! Minimum surface roughness
70   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met = 1)
71
72   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
73   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
74   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn   
75   REAL(wp) ::   rghmin        = -0.28_wp
76   REAL(wp) ::   rgh0          =  0.0329_wp
77   REAL(wp) ::   rghcri        =  0.03_wp
78   REAL(wp) ::   ra1           =  0.92_wp
79   REAL(wp) ::   ra2           =  0.74_wp
80   REAL(wp) ::   rb1           = 16.60_wp
81   REAL(wp) ::   rb2           = 10.10_wp         
82   REAL(wp) ::   re2           =  1.33_wp         
83   REAL(wp) ::   rl1           =  0.107_wp
84   REAL(wp) ::   rl2           =  0.0032_wp
85   REAL(wp) ::   rl3           =  0.0864_wp
86   REAL(wp) ::   rl4           =  0.12_wp
87   REAL(wp) ::   rl5           = 11.9_wp
88   REAL(wp) ::   rl6           =  0.4_wp
89   REAL(wp) ::   rl7           =  0.0_wp
90   REAL(wp) ::   rl8           =  0.48_wp
91   REAL(wp) ::   rm1           =  0.127_wp
92   REAL(wp) ::   rm2           =  0.00336_wp
93   REAL(wp) ::   rm3           =  0.0906_wp
94   REAL(wp) ::   rm4           =  0.101_wp
95   REAL(wp) ::   rm5           = 11.2_wp
96   REAL(wp) ::   rm6           =  0.4_wp
97   REAL(wp) ::   rm7           =  0.0_wp
98   REAL(wp) ::   rm8           =  0.318_wp
99   REAL(wp) ::   rtrans_default =  0.1_wp
100   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
101   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        -
102   REAL(wp) ::   rsbc_psi2, rsbc_psi3, rfact_psi                  !     -           -           -        -
103   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
104   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
105   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        -
106   REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5                     !     -           -           -        -
107   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        -
108   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        -
109
110   !! * Substitutions
111#  include "domzgr_substitute.h90"
112#  include "vectopt_loop_substitute.h90"
113   !!----------------------------------------------------------------------
114   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
115   !! $Id$
116   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
117   !!----------------------------------------------------------------------
118CONTAINS
119
120   INTEGER FUNCTION zdf_gls_alloc()
121      !!----------------------------------------------------------------------
122      !!                ***  FUNCTION zdf_gls_alloc  ***
123      !!----------------------------------------------------------------------
124      ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) ,     &
125         &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     &
126         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     &
127         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), STAT= zdf_gls_alloc )
128         !
129      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
130      IF( zdf_gls_alloc /= 0 )   CALL ctl_warn('zdf_gls_alloc: failed to allocate arrays')
131   END FUNCTION zdf_gls_alloc
132
133
134   SUBROUTINE zdf_gls( kt )
135      !!----------------------------------------------------------------------
136      !!                   ***  ROUTINE zdf_gls  ***
137      !!
138      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
139      !!              coefficients using the GLS turbulent closure scheme.
140      !!----------------------------------------------------------------------
141      INTEGER, INTENT(in) ::   kt ! ocean time step
142      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
143      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars
144      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -
145      REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      -
146      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
147      REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      -
148      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep
149      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar
150      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves
151      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
152      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before
153      REAL(wp), POINTER, DIMENSION(:,:,:) ::   mxlb        ! mixing length at time before
154      REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear
155      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate
156      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
157      REAL(wp), POINTER, DIMENSION(:,:,:) ::   psi         ! psi at time now
158      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a    ! element of the first  matrix diagonal
159      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal
160      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal
161      !!--------------------------------------------------------------------
162      !
163      IF( nn_timing == 1 )  CALL timing_start('zdf_gls')
164      !
165      CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )
166      CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  )
167     
168      ! Preliminary computing
169
170      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
171
172
173      ! variable initialization
174      IF( ln_phioc ) THEN
175         IF( nn_wmix==jp_janssen ) THEN
176            rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking
177            rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking
178            rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking
179         ELSE IF( nn_wmix==jp_craigbanner ) THEN
180            rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf
181            rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf
182            rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp )
183         ENDIF
184      ENDIF
185
186      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
187         avt (:,:,:) = avt_k (:,:,:)
188         avm (:,:,:) = avm_k (:,:,:)
189         avmu(:,:,:) = avmu_k(:,:,:)
190         avmv(:,:,:) = avmv_k(:,:,:) 
191      ENDIF
192
193      ! Compute surface and bottom friction at T-points
194!CDIR NOVERRCHK         
195      DO jj = 2, jpjm1         
196!CDIR NOVERRCHK         
197         DO ji = fs_2, fs_jpim1   ! vector opt.         
198            !
199            ! surface friction
200            ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
201            !   
202            ! bottom friction (explicit before friction)       
203            ! Note that we chose here not to bound the friction as in dynbfr)   
204            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &         
205               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )     
206            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &         
207               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )     
208            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)         
209         END DO         
210      END DO   
211
212      ! Set surface roughness length
213      SELECT CASE ( nn_z0_met )
214      !
215      CASE ( 0 )             ! Constant roughness         
216         zhsro(:,:) = rn_hsro
217      CASE ( 1 )             ! Standard Charnock formula
218         zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro)
219      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
220         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
221         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
222      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
223         WHERE( hsw == 0._wp ) ! surface roughness length according to Charnock formula when sign. wave height 0
224            zhsro = MAX(rn_charn / grav * ustars2, rn_hsro)
225         ELSEWHERE
226            zhsro = MAX(hsw, rn_hsro)
227         END WHERE
228      END SELECT
229
230      ! Compute shear and dissipation rate
231      DO jk = 2, jpkm1
232         DO jj = 2, jpjm1
233            DO ji = fs_2, fs_jpim1   ! vector opt.
234               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
235                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &
236                  &                            / (  fse3uw_n(ji,jj,jk)               &
237                  &                            *    fse3uw_b(ji,jj,jk) )
238               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
239                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
240                  &                            / (  fse3vw_n(ji,jj,jk)               &
241                  &                            *    fse3vw_b(ji,jj,jk) )
242               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
243            END DO
244         END DO
245      END DO
246      !
247      ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
248      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
249
250      ! Save tke at before time step
251      eb  (:,:,:) = en  (:,:,:)
252      mxlb(:,:,:) = mxln(:,:,:)
253
254      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
255         DO jk = 2, jpkm1
256            DO jj = 2, jpjm1 
257               DO ji = fs_2, fs_jpim1   ! vector opt.
258                  zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)
259                  zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbkt(ji,jj)+1) )
260                  zcoef = ( zup / MAX( zdown, rsmall ) )
261                  zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk)
262               END DO
263            END DO
264         END DO
265      ENDIF
266
267      !!---------------------------------!!
268      !!   Equation to prognostic k      !!
269      !!---------------------------------!!
270      !
271      ! Now Turbulent kinetic energy (output in en)
272      ! -------------------------------
273      ! Resolution of a tridiagonal linear system by a "methode de chasse"
274      ! computation from level 2 to jpkm1  (e(1) computed after and e(jpk)=0 ).
275      ! The surface boundary condition are set after
276      ! The bottom boundary condition are also set after. In standard e(bottom)=0.
277      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
278      ! Warning : after this step, en : right hand side of the matrix
279
280      DO jk = 2, jpkm1
281         DO jj = 2, jpjm1
282            DO ji = fs_2, fs_jpim1   ! vector opt.
283               !
284               ! shear prod. at w-point weightened by mask
285               shear(ji,jj,jk) =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
286                  &             + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
287               !
288               ! stratif. destruction
289               buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)
290               !
291               ! shear prod. - stratif. destruction
292               diss = eps(ji,jj,jk)
293               !
294               dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
295               !
296               zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk)          ! production term
297               zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
298               !
299               ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
300               ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style)
301               ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
302               ! Otherwise, this should be rsc_psi/rsc_psi0
303               IF( ln_sigpsi ) THEN
304                  zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
305                  zwall_psi(ji,jj,jk) = rsc_psi /   & 
306                     &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
307               ELSE
308                  zwall_psi(ji,jj,jk) = 1._wp
309               ENDIF
310               !
311               ! building the matrix
312               zcof = rfact_tke * tmask(ji,jj,jk)
313               !
314               ! lower diagonal
315               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
316                  &                      / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
317               !
318               ! upper diagonal
319               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
320                  &                      / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
321               !
322               ! diagonal
323               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
324                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
325               !
326               ! right hand side in en
327               en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
328            END DO
329         END DO
330      END DO
331      !
332      z_elem_b(:,:,jpk) = 1._wp
333      !
334      ! Set surface condition on zwall_psi (1 at the bottom)
335      zwall_psi(:,:,1) = zwall_psi(:,:,2)
336      zwall_psi(:,:,jpk) = 1.
337      !
338      ! Surface boundary condition on tke
339      ! ---------------------------------
340      !
341      SELECT CASE ( nn_bc_surf )
342      !
343      CASE ( 0 )             ! Dirichlet case
344         IF( ln_phioc ) THEN  ! wave induced mixing case with forced/coupled fields
345            IF( nn_wmix==jp_janssen ) THEN
346               ! First level
347               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )
348               z_elem_a(:,:,1) = en(:,:,1)
349               z_elem_c(:,:,1) = 0._wp
350               z_elem_b(:,:,1) = 1._wp
351
352               ! One level below
353               en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
354               z_elem_a(:,:,2) = 0._wp
355               z_elem_c(:,:,2) = 0._wp
356               z_elem_b(:,:,2) = 1._wp
357            ELSE IF( nn_wmix==jp_craigbanner ) THEN
358               en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)
359               en(:,:,1) = MAX(en(:,:,1), rn_emin) 
360               z_elem_a(:,:,1) = en(:,:,1)
361               z_elem_c(:,:,1) = 0._wp
362               z_elem_b(:,:,1) = 1._wp
363               !
364               ! One level below
365               en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) &
366                   &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
367               en(:,:,2) = MAX(en(:,:,2), rn_emin )
368               z_elem_a(:,:,2) = 0._wp 
369               z_elem_c(:,:,2) = 0._wp
370               z_elem_b(:,:,2) = 1._wp
371               !
372            ENDIF
373         ELSE                 ! wave induced mixing case with default values
374            en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp)
375            en(:,:,1) = MAX(en(:,:,1), rn_emin) 
376            z_elem_a(:,:,1) = en(:,:,1)
377            z_elem_c(:,:,1) = 0._wp
378            z_elem_b(:,:,1) = 1._wp
379            !
380            ! One level below
381            en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) &
382                &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
383            en(:,:,2) = MAX(en(:,:,2), rn_emin )
384            z_elem_a(:,:,2) = 0._wp 
385            z_elem_c(:,:,2) = 0._wp
386            z_elem_b(:,:,2) = 1._wp
387            !
388            !
389         ENDIF
390      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
391         IF( ln_phioc ) THEN   ! Shear free case: d(e)/dz=Fw with forced/coupled fields
392            IF( nn_wmix==jp_janssen ) THEN
393               ! Dirichlet conditions at k=1
394               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )
395               z_elem_a(:,:,1) = en(:,:,1)
396               z_elem_c(:,:,1) = 0._wp
397               z_elem_b(:,:,1) = 1._wp
398               ! at k=2, set de/dz=Fw
399               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
400               z_elem_a(:,:,2) = 0._wp
401               zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)
402               en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
403            ELSE IF( nn_wmix==jp_craigbanner ) THEN
404               ! Dirichlet conditions at k=1
405               en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)
406               en(:,:,1)       = MAX(en(:,:,1), rn_emin)     
407               z_elem_a(:,:,1) = en(:,:,1)
408               z_elem_c(:,:,1) = 0._wp
409               z_elem_b(:,:,1) = 1._wp
410               !
411               ! at k=2, set de/dz=Fw
412               !cbr
413               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
414               z_elem_a(:,:,2) = 0._wp
415               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) ))
416               zflxs(:,:)      = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) &
417                    &                           * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
418
419               en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
420               !
421            ENDIF
422         ELSE                  ! Shear free case: d(e)/dz=Fw with default values
423            ! Dirichlet conditions at k=1
424            en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp)
425            en(:,:,1)       = MAX(en(:,:,1), rn_emin)     
426            z_elem_a(:,:,1) = en(:,:,1)
427            z_elem_c(:,:,1) = 0._wp
428            z_elem_b(:,:,1) = 1._wp
429            !
430            ! at k=2, set de/dz=Fw
431            !cbr
432            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
433            z_elem_a(:,:,2) = 0._wp
434            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) ))
435            zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
436                 &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
437
438            en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
439            !
440            !
441         ENDIF
442      END SELECT
443
444      ! Bottom boundary condition on tke
445      ! --------------------------------
446      !
447      SELECT CASE ( nn_bc_bot )
448      !
449      CASE ( 0 )             ! Dirichlet
450         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
451         !                      ! Balance between the production and the dissipation terms
452!CDIR NOVERRCHK
453         DO jj = 2, jpjm1
454!CDIR NOVERRCHK
455            DO ji = fs_2, fs_jpim1   ! vector opt.
456               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
457               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
458               !
459               ! Bottom level Dirichlet condition:
460               z_elem_a(ji,jj,ibot  ) = 0._wp
461               z_elem_c(ji,jj,ibot  ) = 0._wp
462               z_elem_b(ji,jj,ibot  ) = 1._wp
463               en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
464               !
465               ! Just above last level, Dirichlet condition again
466               z_elem_a(ji,jj,ibotm1) = 0._wp
467               z_elem_c(ji,jj,ibotm1) = 0._wp
468               z_elem_b(ji,jj,ibotm1) = 1._wp
469               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
470            END DO
471         END DO
472         !
473      CASE ( 1 )             ! Neumman boundary condition
474         !                     
475!CDIR NOVERRCHK
476         DO jj = 2, jpjm1
477!CDIR NOVERRCHK
478            DO ji = fs_2, fs_jpim1   ! vector opt.
479               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
480               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
481               !
482               ! Bottom level Dirichlet condition:
483               z_elem_a(ji,jj,ibot) = 0._wp
484               z_elem_c(ji,jj,ibot) = 0._wp
485               z_elem_b(ji,jj,ibot) = 1._wp
486               en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
487               !
488               ! Just above last level: Neumann condition
489               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)   ! Remove z_elem_c from z_elem_b
490               z_elem_c(ji,jj,ibotm1) = 0._wp
491            END DO
492         END DO
493         !
494      END SELECT
495
496      ! Matrix inversion (en prescribed at surface and the bottom)
497      ! ----------------------------------------------------------
498      !
499      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
500         DO jj = 2, jpjm1
501            DO ji = fs_2, fs_jpim1    ! vector opt.
502               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
503            END DO
504         END DO
505      END DO
506      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
507         DO jj = 2, jpjm1
508            DO ji = fs_2, fs_jpim1    ! vector opt.
509               z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
510            END DO
511         END DO
512      END DO
513      DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
514         DO jj = 2, jpjm1
515            DO ji = fs_2, fs_jpim1    ! vector opt.
516               en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
517            END DO
518         END DO
519      END DO
520      !                                            ! set the minimum value of tke
521      en(:,:,:) = MAX( en(:,:,:), rn_emin )
522
523      !!----------------------------------------!!
524      !!   Solve prognostic equation for psi    !!
525      !!----------------------------------------!!
526
527      ! Set psi to previous time step value
528      !
529      SELECT CASE ( nn_clos )
530      !
531      CASE( 0 )               ! k-kl  (Mellor-Yamada)
532         DO jk = 2, jpkm1
533            DO jj = 2, jpjm1
534               DO ji = fs_2, fs_jpim1   ! vector opt.
535                  psi(ji,jj,jk)  = eb(ji,jj,jk) * mxlb(ji,jj,jk)
536               END DO
537            END DO
538         END DO
539         !
540      CASE( 1 )               ! k-eps
541         DO jk = 2, jpkm1
542            DO jj = 2, jpjm1
543               DO ji = fs_2, fs_jpim1   ! vector opt.
544                  psi(ji,jj,jk)  = eps(ji,jj,jk)
545               END DO
546            END DO
547         END DO
548         !
549      CASE( 2 )               ! k-w
550         DO jk = 2, jpkm1
551            DO jj = 2, jpjm1
552               DO ji = fs_2, fs_jpim1   ! vector opt.
553                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * mxlb(ji,jj,jk) )
554               END DO
555            END DO
556         END DO
557         !
558      CASE( 3 )               ! generic
559         DO jk = 2, jpkm1
560            DO jj = 2, jpjm1
561               DO ji = fs_2, fs_jpim1   ! vector opt.
562                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 
563               END DO
564            END DO
565         END DO
566         !
567      END SELECT
568      !
569      ! Now gls (output in psi)
570      ! -------------------------------
571      ! Resolution of a tridiagonal linear system by a "methode de chasse"
572      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
573      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
574      ! Warning : after this step, en : right hand side of the matrix
575
576      DO jk = 2, jpkm1
577         DO jj = 2, jpjm1
578            DO ji = fs_2, fs_jpim1   ! vector opt.
579               !
580               ! psi / k
581               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
582               !
583               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
584               dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
585               !
586               rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p
587               !
588               ! shear prod. - stratif. destruction
589               prod = rpsi1 * zratio * shear(ji,jj,jk)
590               !
591               ! stratif. destruction
592               buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )
593               !
594               ! shear prod. - stratif. destruction
595               diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
596               !
597               dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
598               !
599               zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term
600               zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
601               !                                                       
602               ! building the matrix
603               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
604               ! lower diagonal
605               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
606                  &                      / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
607               ! upper diagonal
608               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
609                  &                      / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
610               ! diagonal
611               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
612                  &                       + rdt * zdiss * tmask(ji,jj,jk)
613               !
614               ! right hand side in psi
615               psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
616            END DO
617         END DO
618      END DO
619      !
620      z_elem_b(:,:,jpk) = 1._wp
621
622      ! Surface boundary condition on psi
623      ! ---------------------------------
624      !
625      SELECT CASE ( nn_bc_surf )
626      !
627      CASE ( 0 )             ! Dirichlet boundary conditions
628         IF( ln_phioc ) THEN   ! Wave induced mixing case
629                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
630                               ! balance between the production and the
631                               ! dissipation terms including the wave effect
632            IF( nn_wmix==jp_janssen ) THEN
633               ! Surface value
634               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
635               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
636               z_elem_a(:,:,1) = psi(:,:,1)     
637               z_elem_c(:,:,1) = 0._wp
638               z_elem_b(:,:,1) = 1._wp
639               !
640               ! One level below
641               zex1 = (rmm*ra_sf+rnn)
642               zex2 = (rmm*ra_sf)
643               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2
644               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)
645               z_elem_a(:,:,2) = 0._wp
646               z_elem_c(:,:,2) = 0._wp
647               z_elem_b(:,:,2) = 1._wp
648               !
649               !
650            ELSE IF( nn_wmix==jp_craigbanner ) THEN
651               ! Surface value
652               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
653               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
654               z_elem_a(:,:,1) = psi(:,:,1)
655               z_elem_c(:,:,1) = 0._wp
656               z_elem_b(:,:,1) = 1._wp
657               !
658               ! One level below
659               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) )))
660               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)
661               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
662               z_elem_a(:,:,2) = 0._wp
663               z_elem_c(:,:,2) = 0._wp
664               z_elem_b(:,:,2) = 1._wp
665               !
666               !
667            ENDIF
668         ELSE                  ! Wave induced mixing case with default values
669            ! Surface value
670            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
671            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
672            z_elem_a(:,:,1) = psi(:,:,1)
673            z_elem_c(:,:,1) = 0._wp
674            z_elem_b(:,:,1) = 1._wp
675            !
676            ! One level below
677            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) )))
678            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)
679            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
680            z_elem_a(:,:,2) = 0._wp
681            z_elem_c(:,:,2) = 0._wp
682            z_elem_b(:,:,2) = 1._wp
683            !
684            !
685         ENDIF
686      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
687         IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields
688            IF( nn_wmix==jp_janssen ) THEN
689               !
690               zdep(:,:) = rl_sf * zhsro(:,:)
691               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
692               z_elem_a(:,:,1) = psi(:,:,1)
693               z_elem_c(:,:,1) = 0._wp
694               z_elem_b(:,:,1) = 1._wp
695               !
696               ! Neumann condition at k=2
697               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
698               z_elem_a(:,:,2) = 0._wp
699               !
700               ! Set psi vertical flux at the surface:
701               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)
702               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &
703                  &                   * en(:,:,1)**rmm * zdep
704               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
705               !
706            ELSE IF( nn_wmix==jp_craigbanner ) THEN
707               ! Surface value: Dirichlet
708               zdep(:,:)       = zhsro(:,:) * rl_sf
709               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
710               z_elem_a(:,:,1) = psi(:,:,1)
711               z_elem_c(:,:,1) = 0._wp
712               z_elem_b(:,:,1) = 1._wp
713               !
714               ! Neumann condition at k=2
715               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
716               z_elem_a(:,:,2) = 0._wp
717               !
718               ! Set psi vertical flux at the surface:
719               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
720               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
721               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * &
722                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
723               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
724                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)
725               zflxs(:,:) = zdep(:,:) * zflxs(:,:)
726               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
727               !   
728            ENDIF
729         ELSE                 ! Wave induced mixing case with default values
730            ! Surface value: Dirichlet
731            zdep(:,:)       = zhsro(:,:) * rl_sf
732            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
733            z_elem_a(:,:,1) = psi(:,:,1)
734            z_elem_c(:,:,1) = 0._wp
735            z_elem_b(:,:,1) = 1._wp
736            !
737            ! Neumann condition at k=2
738            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
739            z_elem_a(:,:,2) = 0._wp
740            !
741            ! Set psi vertical flux at the surface:
742            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
743            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
744            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * &
745                       (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
746            zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
747                   & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)
748            zflxs(:,:) = zdep(:,:) * zflxs(:,:)
749            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
750            !   
751            !
752         ENDIF
753      END SELECT
754
755      ! Bottom boundary condition on psi
756      ! --------------------------------
757      !
758      SELECT CASE ( nn_bc_bot )
759      !
760      !
761      CASE ( 0 )             ! Dirichlet
762         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
763         !                      ! Balance between the production and the dissipation terms
764!CDIR NOVERRCHK
765         DO jj = 2, jpjm1
766!CDIR NOVERRCHK
767            DO ji = fs_2, fs_jpim1   ! vector opt.
768               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
769               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
770               zdep(ji,jj) = vkarmn * rn_bfrz0
771               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
772               z_elem_a(ji,jj,ibot) = 0._wp
773               z_elem_c(ji,jj,ibot) = 0._wp
774               z_elem_b(ji,jj,ibot) = 1._wp
775               !
776               ! Just above last level, Dirichlet condition again (GOTM like)
777               zdep(ji,jj) = vkarmn * ( rn_bfrz0 + fse3t(ji,jj,ibotm1) )
778               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
779               z_elem_a(ji,jj,ibotm1) = 0._wp
780               z_elem_c(ji,jj,ibotm1) = 0._wp
781               z_elem_b(ji,jj,ibotm1) = 1._wp
782            END DO
783         END DO
784         !
785      CASE ( 1 )             ! Neumman boundary condition
786         !                     
787!CDIR NOVERRCHK
788         DO jj = 2, jpjm1
789!CDIR NOVERRCHK
790            DO ji = fs_2, fs_jpim1   ! vector opt.
791               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
792               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
793               !
794               ! Bottom level Dirichlet condition:
795               zdep(ji,jj) = vkarmn * rn_bfrz0
796               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
797               !
798               z_elem_a(ji,jj,ibot) = 0._wp
799               z_elem_c(ji,jj,ibot) = 0._wp
800               z_elem_b(ji,jj,ibot) = 1._wp
801               !
802               ! Just above last level: Neumann condition with flux injection
803               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b
804               z_elem_c(ji,jj,ibotm1) = 0.
805               !
806               ! Set psi vertical flux at the bottom:
807               zdep(ji,jj) = rn_bfrz0 + 0.5_wp*fse3t(ji,jj,ibotm1)
808               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   &
809                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
810               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)
811            END DO
812         END DO
813         !
814      END SELECT
815
816      ! Matrix inversion
817      ! ----------------
818      !
819      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
820         DO jj = 2, jpjm1
821            DO ji = fs_2, fs_jpim1    ! vector opt.
822               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
823            END DO
824         END DO
825      END DO
826      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
827         DO jj = 2, jpjm1
828            DO ji = fs_2, fs_jpim1    ! vector opt.
829               z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
830            END DO
831         END DO
832      END DO
833      DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
834         DO jj = 2, jpjm1
835            DO ji = fs_2, fs_jpim1    ! vector opt.
836               psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
837            END DO
838         END DO
839      END DO
840
841      ! Set dissipation
842      !----------------
843
844      SELECT CASE ( nn_clos )
845      !
846      CASE( 0 )               ! k-kl  (Mellor-Yamada)
847         DO jk = 1, jpkm1
848            DO jj = 2, jpjm1
849               DO ji = fs_2, fs_jpim1   ! vector opt.
850                  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)
851               END DO
852            END DO
853         END DO
854         !
855      CASE( 1 )               ! k-eps
856         DO jk = 1, jpkm1
857            DO jj = 2, jpjm1
858               DO ji = fs_2, fs_jpim1   ! vector opt.
859                  eps(ji,jj,jk) = psi(ji,jj,jk)
860               END DO
861            END DO
862         END DO
863         !
864      CASE( 2 )               ! k-w
865         DO jk = 1, jpkm1
866            DO jj = 2, jpjm1
867               DO ji = fs_2, fs_jpim1   ! vector opt.
868                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
869               END DO
870            END DO
871         END DO
872         !
873      CASE( 3 )               ! generic
874         zcoef = rc0**( 3._wp  + rpp/rnn )
875         zex1  =      ( 1.5_wp + rmm/rnn )
876         zex2  = -1._wp / rnn
877         DO jk = 1, jpkm1
878            DO jj = 2, jpjm1
879               DO ji = fs_2, fs_jpim1   ! vector opt.
880                  eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
881               END DO
882            END DO
883         END DO
884         !
885      END SELECT
886
887      ! Limit dissipation rate under stable stratification
888      ! --------------------------------------------------
889      DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time
890         DO jj = 2, jpjm1
891            DO ji = fs_2, fs_jpim1    ! vector opt.
892               ! limitation
893               eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
894               mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
895               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
896               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
897               IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
898            END DO
899         END DO
900      END DO 
901
902      !
903      ! Stability function and vertical viscosity and diffusivity
904      ! ---------------------------------------------------------
905      !
906      SELECT CASE ( nn_stab_func )
907      !
908      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
909         DO jk = 2, jpkm1
910            DO jj = 2, jpjm1
911               DO ji = fs_2, fs_jpim1   ! vector opt.
912                  ! zcof =  l²/q²
913                  zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
914                  ! Gh = -N²l²/q²
915                  gh = - rn2(ji,jj,jk) * zcof
916                  gh = MIN( gh, rgh0   )
917                  gh = MAX( gh, rghmin )
918                  ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
919                  sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) )
920                  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)
921                  !
922                  ! Store stability function in avmu and avmv
923                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
924                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
925               END DO
926            END DO
927         END DO
928         !
929      CASE ( 2, 3 )               ! Canuto stability functions
930         DO jk = 2, jpkm1
931            DO jj = 2, jpjm1
932               DO ji = fs_2, fs_jpim1   ! vector opt.
933                  ! zcof =  l²/q²
934                  zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
935                  ! Gh = -N²l²/q²
936                  gh = - rn2(ji,jj,jk) * zcof
937                  gh = MIN( gh, rgh0   )
938                  gh = MAX( gh, rghmin )
939                  gh = gh * rf6
940                  ! Gm =  M²l²/q² Shear number
941                  shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall )
942                  gm = MAX( shr * zcof , 1.e-10 )
943                  gm = gm * rf6
944                  gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm )
945                  ! Stability functions from Canuto
946                  rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm
947                  sm = (rs0 - rs1*gh + rs2*gm) / rcff
948                  sh = (rs4 - rs5*gh + rs6*gm) / rcff
949                  !
950                  ! Store stability function in avmu and avmv
951                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
952                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
953               END DO
954            END DO
955         END DO
956         !
957      END SELECT
958
959      ! Boundary conditions on stability functions for momentum (Neumann):
960      ! Lines below are useless if GOTM style Dirichlet conditions are used
961
962      avmv(:,:,1) = avmv(:,:,2)
963
964      DO jj = 2, jpjm1
965         DO ji = fs_2, fs_jpim1   ! vector opt.
966            avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj))
967         END DO
968      END DO
969
970      ! Compute diffusivities/viscosities
971      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
972      DO jk = 1, jpk
973         DO jj = 2, jpjm1
974            DO ji = fs_2, fs_jpim1   ! vector opt.
975               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)
976               zav           = zsqen * avmu(ji,jj,jk)
977               avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine
978               zav           = zsqen * avmv(ji,jj,jk)
979               avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom
980            END DO
981         END DO
982      END DO
983      !
984      ! Lateral boundary conditions (sign unchanged)
985      avt(:,:,1)  = 0._wp
986      CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
987
988      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
989         DO jj = 2, jpjm1
990            DO ji = fs_2, fs_jpim1   ! vector opt.
991               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
992               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
993            END DO
994         END DO
995      END DO
996      avmu(:,:,1) = 0._wp             ;   avmv(:,:,1) = 0._wp                 ! set surface to zero
997      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )       ! Lateral boundary conditions
998
999      IF(ln_ctl) THEN
1000         CALL prt_ctl( tab3d_1=en  , clinfo1=' gls  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
1001         CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls  - u: ', mask1=umask,                   &
1002            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
1003      ENDIF
1004      !
1005      avt_k (:,:,:) = avt (:,:,:)
1006      avm_k (:,:,:) = avm (:,:,:)
1007      avmu_k(:,:,:) = avmu(:,:,:)
1008      avmv_k(:,:,:) = avmv(:,:,:)
1009      !
1010      CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro )
1011      CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
1012      !
1013      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls')
1014      !
1015      !
1016   END SUBROUTINE zdf_gls
1017
1018
1019   SUBROUTINE zdf_gls_init
1020      !!----------------------------------------------------------------------
1021      !!                  ***  ROUTINE zdf_gls_init  ***
1022      !!                     
1023      !! ** Purpose :   Initialization of the vertical eddy diffivity and
1024      !!      viscosity when using a gls turbulent closure scheme
1025      !!
1026      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
1027      !!      called at the first timestep (nit000)
1028      !!
1029      !! ** input   :   Namlist namzdf_gls
1030      !!
1031      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
1032      !!
1033      !!----------------------------------------------------------------------
1034      USE dynzdf_exp
1035      USE trazdf_exp
1036      !
1037      INTEGER ::   jk    ! dummy loop indices
1038      INTEGER ::   ios   ! Local integer output status for namelist read
1039      REAL(wp)::   zcr   ! local scalar
1040      !!
1041      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
1042         &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
1043         &            rn_crban_default, rn_charn, rn_frac_hs,&
1044         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
1045         &            nn_stab_func, nn_clos
1046      !!----------------------------------------------------------
1047      !
1048      IF( nn_timing == 1 )  CALL timing_start('zdf_gls_init')
1049      !
1050      REWIND( numnam_ref )              ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
1051      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
1052901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in reference namelist', lwp )
1053
1054      REWIND( numnam_cfg )              ! Namelist namzdf_gls in configuration namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
1055      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
1056902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp )
1057      IF(lwm) WRITE ( numond, namzdf_gls )
1058
1059      IF(lwp) THEN                     !* Control print
1060         WRITE(numout,*)
1061         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
1062         WRITE(numout,*) '~~~~~~~~~~~~'
1063         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
1064         WRITE(numout,*) '      minimum value of en                           rn_emin        = ', rn_emin
1065         WRITE(numout,*) '      minimum value of eps                          rn_epsmin      = ', rn_epsmin
1066         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim  = ', ln_length_lim
1067         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp   = ', rn_clim_galp
1068         WRITE(numout,*) '      TKE Surface boundary condition                nn_bc_surf     = ', nn_bc_surf
1069         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
1070         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
1071         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default
1072         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
1073         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
1074         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs
1075         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
1076         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
1077         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
1078         WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0
1079      ENDIF
1080
1081      !                                !* allocate gls arrays
1082      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
1083
1084      !                                !* Check of some namelist values
1085      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' ) 
1086      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' ) 
1087      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' ) 
1088      IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' ) 
1089      IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN
1090         CALL ctl_warn('W A R N I N G:  ln_rough=.TRUE. but nn_z0_met is not 3 - resetting nn_z0_met to 3')
1091         nn_z0_met = 3
1092      ENDIF
1093      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' ) 
1094      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
1095
1096      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
1097      !
1098      CASE( 0 )                              ! k-kl  (Mellor-Yamada)
1099         !
1100         IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
1101         rpp     = 0._wp
1102         rmm     = 1._wp
1103         rnn     = 1._wp
1104         rsc_tke = 1.96_wp
1105         rsc_psi = 1.96_wp
1106         rpsi1   = 0.9_wp
1107         rpsi3p  = 1._wp
1108         rpsi2   = 0.5_wp
1109         !
1110         SELECT CASE ( nn_stab_func )
1111         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
1112         CASE( 2 )      ;   rpsi3m = 2.62_wp       ! Canuto A stability functions
1113         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
1114         END SELECT
1115         !
1116      CASE( 1 )                              ! k-eps
1117         !
1118         IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
1119         rpp     =  3._wp
1120         rmm     =  1.5_wp
1121         rnn     = -1._wp
1122         rsc_tke =  1._wp
1123         rsc_psi =  1.2_wp  ! Schmidt number for psi
1124         rpsi1   =  1.44_wp
1125         rpsi3p  =  1._wp
1126         rpsi2   =  1.92_wp
1127         !
1128         SELECT CASE ( nn_stab_func )
1129         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
1130         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
1131         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
1132         END SELECT
1133         !
1134      CASE( 2 )                              ! k-omega
1135         !
1136         IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
1137         rpp     = -1._wp
1138         rmm     =  0.5_wp
1139         rnn     = -1._wp
1140         rsc_tke =  2._wp
1141         rsc_psi =  2._wp
1142         rpsi1   =  0.555_wp
1143         rpsi3p  =  1._wp
1144         rpsi2   =  0.833_wp
1145         !
1146         SELECT CASE ( nn_stab_func )
1147         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
1148         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
1149         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
1150         END SELECT
1151         !
1152      CASE( 3 )                              ! generic
1153         !
1154         IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
1155         rpp     = 2._wp
1156         rmm     = 1._wp
1157         rnn     = -0.67_wp
1158         rsc_tke = 0.8_wp
1159         rsc_psi = 1.07_wp
1160         rpsi1   = 1._wp
1161         rpsi3p  = 1._wp
1162         rpsi2   = 1.22_wp
1163         !
1164         SELECT CASE ( nn_stab_func )
1165         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
1166         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
1167         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
1168         END SELECT
1169         !
1170      END SELECT
1171
1172      !
1173      SELECT CASE ( nn_stab_func )     !* set the parameters of the stability functions
1174      !
1175      CASE ( 0 )                             ! Galperin stability functions
1176         !
1177         IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
1178         rc2     =  0._wp
1179         rc3     =  0._wp
1180         rc_diff =  1._wp
1181         rc0     =  0.5544_wp
1182         rcm_sf  =  0.9884_wp
1183         rghmin  = -0.28_wp
1184         rgh0    =  0.0233_wp
1185         rghcri  =  0.02_wp
1186         !
1187      CASE ( 1 )                             ! Kantha-Clayson stability functions
1188         !
1189         IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
1190         rc2     =  0.7_wp
1191         rc3     =  0.2_wp
1192         rc_diff =  1._wp
1193         rc0     =  0.5544_wp
1194         rcm_sf  =  0.9884_wp
1195         rghmin  = -0.28_wp
1196         rgh0    =  0.0233_wp
1197         rghcri  =  0.02_wp
1198         !
1199      CASE ( 2 )                             ! Canuto A stability functions
1200         !
1201         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
1202         rs0 = 1.5_wp * rl1 * rl5*rl5
1203         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
1204         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
1205         rs4 = 2._wp * rl5
1206         rs5 = 2._wp * rl4
1207         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
1208            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
1209         rd0 = 3._wp * rl5*rl5
1210         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
1211         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
1212         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
1213         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
1214         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
1215         rc0 = 0.5268_wp
1216         rf6 = 8._wp / (rc0**6._wp)
1217         rc_diff = SQRT(2._wp) / (rc0**3._wp)
1218         rcm_sf  =  0.7310_wp
1219         rghmin  = -0.28_wp
1220         rgh0    =  0.0329_wp
1221         rghcri  =  0.03_wp
1222         !
1223      CASE ( 3 )                             ! Canuto B stability functions
1224         !
1225         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
1226         rs0 = 1.5_wp * rm1 * rm5*rm5
1227         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
1228         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
1229         rs4 = 2._wp * rm5
1230         rs5 = 2._wp * rm4
1231         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)
1232         rd0 = 3._wp * rm5*rm5
1233         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
1234         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
1235         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
1236         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
1237         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
1238         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
1239         rf6 = 8._wp / ( rc0**6._wp )
1240         rc_diff = SQRT(2._wp)/(rc0**3.)
1241         rcm_sf  =  0.7470_wp
1242         rghmin  = -0.28_wp
1243         rgh0    =  0.0444_wp
1244         rghcri  =  0.0414_wp
1245         !
1246      END SELECT
1247   
1248      !                                !* Set Schmidt number for psi diffusion in the wave breaking case
1249      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1250      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1251      IF( ln_sigpsi ) THEN
1252         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
1253         ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
1254         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
1255         ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn
1256         rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.))
1257      ELSE
1258         rsc_psi0 = rsc_psi
1259      ENDIF
1260 
1261      !                                !* Shear free turbulence parameters
1262      !
1263      ra_sf  = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) &
1264               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
1265
1266      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN
1267         rl_sf = vkarmn
1268      ELSE
1269         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke          &
1270                 &                                       + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) &
1271                 &                                                *SQRT(rsc_tke*(rsc_tke                 &
1272                 &                                                   + 24._wp*rsc_psi0*rpsi2)) )         &
1273                 &                                         /(12._wp*rnn**2.)                             &
1274                 &                                       )
1275      ENDIF
1276
1277      !
1278      IF(lwp) THEN                     !* Control print
1279         WRITE(numout,*)
1280         WRITE(numout,*) 'Limit values'
1281         WRITE(numout,*) '~~~~~~~~~~~~'
1282         WRITE(numout,*) 'Parameter  m = ',rmm
1283         WRITE(numout,*) 'Parameter  n = ',rnn
1284         WRITE(numout,*) 'Parameter  p = ',rpp
1285         WRITE(numout,*) 'rpsi1   = ',rpsi1
1286         WRITE(numout,*) 'rpsi2   = ',rpsi2
1287         WRITE(numout,*) 'rpsi3m  = ',rpsi3m
1288         WRITE(numout,*) 'rpsi3p  = ',rpsi3p
1289         WRITE(numout,*) 'rsc_tke = ',rsc_tke
1290         WRITE(numout,*) 'rsc_psi = ',rsc_psi
1291         WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0
1292         WRITE(numout,*) 'rc0     = ',rc0
1293         WRITE(numout,*)
1294         WRITE(numout,*) 'Shear free turbulence parameters:'
1295         WRITE(numout,*) 'rcm_sf  = ',rcm_sf
1296         WRITE(numout,*) 'ra_sf   = ',ra_sf
1297         WRITE(numout,*) 'rl_sf   = ',rl_sf
1298         WRITE(numout,*)
1299      ENDIF
1300
1301      !                                !* Constants initialization
1302      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
1303      rc03  = rc02 * rc0
1304      rc04  = rc03 * rc0
1305      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf
1306      rsbc_tke2         = rdt * rn_crban_default / rl_sf
1307      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp )
1308      rtrans_default    = 0.2_wp / zcr
1309      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
1310      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
1311      rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
1312      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
1313
1314      rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
1315      rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke
1316
1317      !                                !* Wall proximity function
1318      zwall (:,:,:) = 1._wp * tmask(:,:,:)
1319
1320      !                                !* set vertical eddy coef. to the background value
1321      DO jk = 1, jpk
1322         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1323         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
1324         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1325         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1326      END DO
1327      !                             
1328      CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
1329      !
1330      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init')
1331      !
1332   END SUBROUTINE zdf_gls_init
1333
1334
1335   SUBROUTINE gls_rst( kt, cdrw )
1336      !!---------------------------------------------------------------------
1337      !!                   ***  ROUTINE ts_rst  ***
1338      !!                     
1339      !! ** Purpose :   Read or write TKE file (en) in restart file
1340      !!
1341      !! ** Method  :   use of IOM library
1342      !!                if the restart does not contain TKE, en is either
1343      !!                set to rn_emin or recomputed (nn_igls/=0)
1344      !!----------------------------------------------------------------------
1345      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1346      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1347      !
1348      INTEGER ::   jit, jk   ! dummy loop indices
1349      INTEGER ::   id1, id2, id3, id4, id5, id6
1350      INTEGER ::   ji, jj, ikbu, ikbv
1351      REAL(wp)::   cbx, cby
1352      !!----------------------------------------------------------------------
1353      !
1354      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1355         !                                   ! ---------------
1356         IF( ln_rstart ) THEN                   !* Read the restart file
1357            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
1358            id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
1359            id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
1360            id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
1361            id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
1362            id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
1363            !
1364            IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
1365               CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
1366               CALL iom_get( numror, jpdom_autoglo, 'avt'   , avt    )
1367               CALL iom_get( numror, jpdom_autoglo, 'avm'   , avm    )
1368               CALL iom_get( numror, jpdom_autoglo, 'avmu'  , avmu   )
1369               CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   )
1370               CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
1371            ELSE                       
1372               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
1373               en  (:,:,:) = rn_emin
1374               mxln(:,:,:) = 0.05       
1375               avt_k (:,:,:) = avt (:,:,:)
1376               avm_k (:,:,:) = avm (:,:,:)
1377               avmu_k(:,:,:) = avmu(:,:,:)
1378               avmv_k(:,:,:) = avmv(:,:,:)
1379               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
1380            ENDIF
1381         ELSE                                   !* Start from rest
1382            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
1383            en  (:,:,:) = rn_emin
1384            mxln(:,:,:) = 0.05       
1385         ENDIF
1386         !
1387      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1388         !                                   ! -------------------
1389         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1390         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
1391         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
1392         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
1393         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
1394         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
1395         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   )
1396         !
1397      ENDIF
1398      !
1399   END SUBROUTINE gls_rst
1400
1401#else
1402   !!----------------------------------------------------------------------
1403   !!   Dummy module :                                        NO TKE scheme
1404   !!----------------------------------------------------------------------
1405   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfgls = .FALSE.   !: TKE flag
1406CONTAINS
1407   SUBROUTINE zdf_gls_init           ! Empty routine
1408      WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
1409   END SUBROUTINE zdf_gls_init
1410   SUBROUTINE zdf_gls( kt )          ! Empty routine
1411      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
1412   END SUBROUTINE zdf_gls
1413   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
1414      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1415      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1416      WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw
1417   END SUBROUTINE gls_rst
1418#endif
1419
1420   !!======================================================================
1421END MODULE zdfgls
1422
Note: See TracBrowser for help on using the repository browser.