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/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 2299

Last change on this file since 2299 was 2299, checked in by cbricaud, 14 years ago

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