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

Last change on this file since 2395 was 2395, checked in by rblod, 13 years ago

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