source: branches/DEV_r1784_GLS/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 2048

Last change on this file since 2048 was 2048, checked in by cbricaud, 11 years ago

add GLS Vertical scheme

File size: 57.5 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) ::   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   REAL(wp) ::   rn_cm_sf   = 0.73_wp         ! Shear free turbulence parameters
62   REAL(wp) ::   rn_a_sf    = -2.0_wp         ! Must be negative -2 < rn_a_sf < -1
63   REAL(wp) ::   rn_l_sf    =  0.2_wp         ! 0 <rn_l_sf<vkarmn   
64   REAL(wp) ::   rn_ghmin   = -0.28
65   REAL(wp) ::   rn_gh0     = 0.0329
66   REAL(wp) ::   rn_ghcri   = 0.03
67   REAL(wp) ::   rn_a1      =  0.92_wp
68   REAL(wp) ::   rn_a2      =  0.74_wp
69   REAL(wp) ::   rn_b1      = 16.60_wp
70   REAL(wp) ::   rn_b2      = 10.10_wp         
71   REAL(wp) ::   rn_e1      =  1.80_wp         
72   REAL(wp) ::   rn_e2      =  1.33_wp         
73   REAL(wp) ::   rn_l1      =  0.107_wp
74   REAL(wp) ::   rn_l2      =  0.0032_wp
75   REAL(wp) ::   rn_l3      =  0.0864_wp
76   REAL(wp) ::   rn_l4      =  0.12_wp
77   REAL(wp) ::   rn_l5      = 11.9_wp
78   REAL(wp) ::   rn_l6      =  0.4_wp
79   REAL(wp) ::   rn_l7      =  0.0_wp
80   REAL(wp) ::   rn_l8      =  0.48_wp
81   REAL(wp) ::   rn_m1      =  0.127_wp
82   REAL(wp) ::   rn_m2      =  0.00336_wp
83   REAL(wp) ::   rn_m3      =  0.0906_wp
84   REAL(wp) ::   rn_m4      =  0.101_wp
85   REAL(wp) ::   rn_m5      = 11.2_wp
86   REAL(wp) ::   rn_m6      =  0.4_wp
87   REAL(wp) ::   rn_m7      =  0.0_wp
88   REAL(wp) ::   rn_m8      =  0.318_wp
89   REAL(wp) ::   rn_c02
90   REAL(wp) ::   rn_c02r
91   REAL(wp) ::   rn_c03
92   REAL(wp) ::   rn_c04
93   REAL(wp) ::   rn_c03_sqrt2_galp
94   REAL(wp) ::   rn_sbc_mb
95   REAL(wp) ::   rn_sbc_std
96   REAL(wp) ::   rn_sbc_tke1
97   REAL(wp) ::   rn_sbc_tke2
98   REAL(wp) ::   rn_sbc_tke3
99   REAL(wp) ::   rn_sbc_psi1
100   REAL(wp) ::   rn_sbc_psi2
101   REAL(wp) ::   rn_sbc_psi3
102   REAL(wp) ::   rn_sbc_zs
103   REAL(wp) ::   zfact_tke, zfact_psi
104   REAL(wp) ::   rn_c0, rn_c2, rn_c3, rn_f6, rn_cff, rn_c_diff
105   REAL(wp) ::   rn_s0, rn_s1, rn_s2, rn_s4, rn_s5, rn_s6
106   REAL(wp) ::   rn_d0, rn_d1, rn_d2, rn_d3, rn_d4, rn_d5
107   REAL(wp) ::   rn_sc_tke, rn_sc_psi, rn_psi1, rn_psi2, rn_psi3, rn_sc_psi0
108   REAL(wp) ::   rn_psi3m, rn_psi3p, zpp, zmm, znn, epslim
109
110   !! * Substitutions
111#  include "domzgr_substitute.h90"
112#  include "vectopt_loop_substitute.h90"
113   !!----------------------------------------------------------------------
114   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
115   !! $Id: zdfgls.F90 1201 2008-09-24 13:24:21Z rblod $
116   !! Software governed by the CeCILL licence (modipsl/doc/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) = rn_sbc_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(rn_sbc_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)  = rn_c03 * 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. + rn_e2 * ( 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 rn_sc_psi*zwall/rn_sc_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 rn_sc_psi/rn_sc_psi0
262               IF (ln_sigpsi) THEN
263                 zsigpsi = MIN(1., zesh2/eps(ji,jj,jk)) ! 0. <= zsigpsi <= 1.
264                 zwall_psi(ji,jj,jk) = rn_sc_psi / (zsigpsi * rn_sc_psi + (1.-zsigpsi) * rn_sc_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 = zfact_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) = rn_sc_psi/rn_sc_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( rn_sbc_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( rn_sbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+gdepw(:,:,2))/zhsro(:,:) )**rn_a_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( rn_c02r * 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( rn_c02r * 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( rn_sbc_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(:,:) = rn_sbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+gdept(:,:,1))/zhsro(:,:) )**(1.5*rn_a_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( rn_c02r * 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( rn_c02r * 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( rn_c02r * 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( rn_c02r * 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)) / (rn_c0 * 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)  = rn_c02 * en(ji,jj,jk) * mxln(ji,jj,jk)**znn 
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               rn_psi3 = dir*rn_psi3m+(1-dir)*rn_psi3p
519               !
520               ! shear prod. - stratif. destruction
521               prod = rn_psi1 * zratio * shear(ji,jj,jk)
522               !
523               ! stratif. destruction
524               buoy = rn_psi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk))
525               !
526               ! shear prod. - stratif. destruction
527               diss = rn_psi2 * 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 = zfact_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(:,:) = rn_l_sf * zhsro(:,:)
566      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * 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(:,:) + gdepw(:,:,2))**(zmm*rn_a_sf+znn) )  &
573        &                       / zhsro(:,:)**(zmm*rn_a_sf)
574      psi (:,:,2) = rn_sbc_psi1 * ustars2(:,:)**zmm * 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) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * 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(:,:) + gdepw(:,:,2) )
591      psi (:,:,2) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * 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(:,:) = rn_l_sf * zhsro(:,:)
603      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * 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(:,:) + gdept(:,:,1))**(zmm*rn_a_sf+znn-1.) / zhsro(:,:)**(zmm*rn_a_sf)
614      zflxs(:,:) = rn_sbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
615                 & * en(:,:,1)**zmm * zdep         
616      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
617      !
618      ELSE                   ! No wave induced mixing
619      !
620      zdep(:,:) = vkarmn * zhsro(:,:)
621      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * 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(:,:) + gdept(:,:,1)
632      zflxs(:,:) = rn_sbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**zmm * zdep**(znn-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) = rn_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn
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) = rn_c0**zpp * en(ji,jj,ibot  )**zmm * zdep(ji,jj)**znn
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) = rn_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn
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 = rn_sbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * &
694              & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**zmm * zdep(ji,jj)**(znn-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) = rn_c03 * 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) = rn_c04 * 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) = rn_c0**(3.+zpp/znn) * en(ji,jj,jk)**(1.5+zmm/znn) * psi(ji,jj,jk)**(-1./znn)
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)  = rn_c03 * 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( 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, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri))
804               gh = MIN( gh, rn_gh0   )
805               gh = MAX( gh, rn_ghmin )
806               ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
807               sh = rn_a2*( 1.-6.*rn_a1/rn_b1 ) / ( 1.-3.*rn_a2*gh*(6.*rn_a1+rn_b2*( 1.-rn_c3 ) ) )
808               sm = ( rn_b1**(-1./3.) + ( 18*rn_a1*rn_a1 + 9.*rn_a1*rn_a2*(1.-rn_c2) )*sh*gh ) / (1.-9.*rn_a1*rn_a2*gh)
809               !
810               ! Store stability function in avmu and avmv
811               avmu(ji,jj,jk) = rn_c_diff * sh * tmask(ji,jj,jk)
812               avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk)
813            END DO
814         END DO
815      END DO
816      !
817      CASE ( 2, 3 )               ! Canuto stability functions
818      !
819      DO jk = 2, jpkm1
820         DO jj = 2, jpjm1
821            DO ji = fs_2, fs_jpim1   ! vector opt.
822               ! zcof =  l²/q²
823               zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) )
824               ! Gh = -N²l²/q²
825               gh = - rn2(ji,jj,jk) * zcof
826!               gh = MIN(gh, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri))
827               gh = MIN( gh, rn_gh0   )
828               gh = MAX( gh, rn_ghmin )
829               gh = gh*rn_f6
830               ! Gm =  M²l²/q² Shear number
831               shr = shear(ji,jj,jk)/MAX(avm(ji,jj,jk), rsmall)
832               gm = shr * zcof
833               gm = MAX(gm, 1.e-10)
834               gm = gm*rn_f6
835               gm = MIN ( (rn_d0 - rn_d1*gh + rn_d3*gh**2.)/(rn_d2-rn_d4*gh) , gm )
836               ! Stability functions from Canuto
837               rn_cff = rn_d0 - rn_d1*gh +rn_d2*gm + rn_d3*gh**2. - rn_d4*gh*gm + rn_d5*gm**2.
838               sm = (rn_s0 - rn_s1*gh + rn_s2*gm) / rn_cff
839               sh = (rn_s4 - rn_s5*gh + rn_s6*gm) / rn_cff
840               !
841               ! Store stability function in avmu and avmv
842               avmu(ji,jj,jk) = rn_c_diff * sh * tmask(ji,jj,jk)
843               avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk)
844            END DO
845         END DO
846      END DO
847      !
848      END SELECT
849
850      ! Boundary conditions on stability functions for momentum (Neumann):
851      ! Lines below are useless if GOTM style Dirichlet conditions are used
852      DO jj = 2, jpjm1
853         DO ji = fs_2, fs_jpim1   ! vector opt.
854            avmv(ji,jj,1) = rn_cm_sf / SQRT(2.)
855         END DO
856      END DO
857      DO jj = 2, jpjm1
858         DO ji = fs_2, fs_jpim1   ! vector opt.
859            ibot=mbathy(ji,jj)
860            ibotm1=ibot-1
861            avmv(ji,jj,ibot) = rn_c0 / SQRT(2.)
862         END DO
863      END DO
864
865      ! Compute diffusivities/viscosities
866      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
867      DO jk = 1, jpk
868         DO jj = 2, jpjm1
869            DO ji = fs_2, fs_jpim1   ! vector opt.
870               zsqen           = SQRT(2.*en(ji,jj,jk)) * mxln(ji,jj,jk)
871               zav             = zsqen * avmu(ji,jj,jk)
872               avt(ji,jj,jk)   = MAX(zav,avtb(jk))*tmask(ji,jj,jk) ! apply mask for zdfmxl routine
873               zav             = zsqen * avmv(ji,jj,jk)
874               avm(ji,jj,jk)   = MAX(zav,avmb(jk)) ! Note that avm is not masked at the surface and the bottom
875            END DO
876         END DO
877      END DO
878
879      !
880      ! Lateral boundary conditions (sign unchanged)
881      !
882      avt(:,:,1)  = 0.
883      CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
884
885      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
886         DO jj = 2, jpjm1
887            DO ji = fs_2, fs_jpim1   ! vector opt.
888               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)   &
889               &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
890               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)   &
891               &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
892            END DO
893         END DO
894      END DO
895
896      avmu(:,:,1) = 0.
897      avmv(:,:,1) = 0.
898
899      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
900
901      IF(ln_ctl) THEN
902         CALL prt_ctl( tab3d_1=en  , clinfo1=' gls  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
903         CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls  - u: ', mask1=umask,                   &
904            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
905      ENDIF
906      !
907   END SUBROUTINE zdf_gls
908
909   SUBROUTINE zdf_gls_init
910      !!----------------------------------------------------------------------
911      !!                  ***  ROUTINE zdf_gls_init  ***
912      !!                     
913      !! ** Purpose :   Initialization of the vertical eddy diffivity and
914      !!      viscosity when using a gls turbulent closure scheme
915      !!
916      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
917      !!      called at the first timestep (nit000)
918      !!
919      !! ** input   :   Namlist namzdf_gls
920      !!
921      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
922      !!
923      !!----------------------------------------------------------------------
924      USE dynzdf_exp
925      USE trazdf_exp
926      !
927# if defined key_vectopt_memory
928      INTEGER ::   ji, jj, jk   ! dummy loop indices
929# else
930      INTEGER ::           jk   ! dummy loop indices
931# endif
932      REAL(wp)::   zcr
933      !!
934      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
935         &            clim_galp, ln_crban, ln_sigpsi,    &
936         &            rn_crban, rn_charn,                &
937         &            nn_tkebc_surf, nn_tkebc_bot,       &
938         &            nn_psibc_surf, nn_psibc_bot,       &
939         &            nn_stab_func, nn_clos
940      !!----------------------------------------------------------
941
942      ! Read Namelist namzdf_gls
943      ! ------------------------
944      REWIND ( numnam )
945      READ   ( numnam, namzdf_gls )
946
947      ! Parameter control and print
948      ! ---------------------------
949      ! Control print
950      IF(lwp) THEN
951         WRITE(numout,*)
952         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
953         WRITE(numout,*) '~~~~~~~~~~~~'
954         WRITE(numout,*) '          Namelist namzdf_gls : set gls mixing parameters'
955         WRITE(numout,*) '             minimum value of en                       rn_emin  = ', rn_emin
956         WRITE(numout,*) '             minimum value of eps                    rn_epsmin  = ', rn_epsmin
957         WRITE(numout,*) '             Surface roughness (m)                         hsro = ', hsro
958         WRITE(numout,*) '             Bottom roughness (m)                        rn_hbro = ', rn_hbro
959         WRITE(numout,*) '             Limit dissipation rate under stable stratification ln_length_lim = ',ln_length_lim
960         WRITE(numout,*) '             Galperin length scale limitation coef (Standard: 0.53, Holt: 0.26) clim_galp = ', clim_galp
961         WRITE(numout,*) '             TKE Surface boundary condition       nn_tkebc_surf = ', nn_tkebc_surf
962         WRITE(numout,*) '             TKE Bottom boundary condition        nn_tkebc_bot  = ', nn_tkebc_bot
963         WRITE(numout,*) '             PSI Surface boundary condition       nn_psibc_surf = ', nn_psibc_surf
964         WRITE(numout,*) '             PSI Bottom boundary condition        nn_psibc_bot  = ', nn_psibc_bot
965         WRITE(numout,*) '             Craig and Banner scheme                  ln_crban  = ', ln_crban
966         WRITE(numout,*) '             Modify psi Schmidt number (wb case)     ln_sigpsi  = ', ln_sigpsi
967         WRITE(numout,*) '             Craig and Banner coef.                   rn_crban  = ', rn_crban
968         WRITE(numout,*) '             Charnock coef.                           rn_charn  = ', rn_charn
969         WRITE(numout,*) '             Stability functions                   nn_stab_func = ',nn_stab_func
970         WRITE(numout,*) '             Closure                                 nn_clos    = ',nn_clos
971         WRITE(numout,*)
972      ENDIF
973
974      ! Check of some namelist values
975      IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' )
976      IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' )
977      IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' )
978      IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' )
979      IF( nn_stab_func < 0 .OR. nn_stab_func > 3) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
980      IF( nn_clos < 0 .OR. nn_clos > 3) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
981
982      ! Initialisation of the parameters for the choosen closure
983      ! --------------------------------------------------------
984      !
985      SELECT CASE ( nn_clos )
986      !
987      CASE( 0 )               ! k-kl  (Mellor-Yamada)
988      !
989      IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
990      zpp       = 0.
991      zmm       = 1.
992      znn       = 1.
993      rn_sc_tke = 1.96
994      rn_sc_psi = 1.96
995      rn_psi1 = 0.9
996      rn_psi3p = 1.
997      rn_psi2  = 0.5
998      !
999         SELECT CASE ( nn_stab_func )
1000         !
1001         CASE( 0, 1 )               ! G88 or KC stability functions
1002            rn_psi3m = 2.53
1003         CASE( 2 )                  ! Canuto A stability functions
1004            rn_psi3m = 2.38
1005         CASE( 3 )                  ! Canuto B stability functions
1006            rn_psi3m = 2.38         ! caution : constant not identified
1007         END SELECT
1008      !
1009      CASE( 1 )               ! k-eps
1010      !
1011      IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
1012      zpp       = 3.
1013      zmm       = 1.5
1014      znn       = -1.
1015      rn_sc_tke = 1.
1016      rn_sc_psi = 1.3  ! Schmidt number for psi
1017      rn_psi1 = 1.44
1018      rn_psi3p = 1.
1019      rn_psi2  = 1.92
1020      !
1021        SELECT CASE ( nn_stab_func )
1022        !
1023        CASE( 0, 1 )               ! G88 or KC stability functions
1024           rn_psi3m = -0.52
1025        CASE( 2 )                  ! Canuto A stability functions
1026           rn_psi3m = -0.629 
1027        CASE( 3 )                  ! Canuto B stability functions
1028           rn_psi3m = -0.566
1029        END SELECT
1030      !
1031      CASE( 2 )               ! k-omega
1032      !
1033      IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
1034      zpp       = -1.
1035      zmm       = 0.5
1036      znn       = -1.
1037      rn_sc_tke = 2.
1038      rn_sc_psi = 2.
1039      rn_psi1 = 0.555
1040      rn_psi3p = 1.
1041      rn_psi2  = 0.833
1042      !
1043        SELECT CASE ( nn_stab_func )
1044        !
1045        CASE( 0, 1 )               ! G88 or KC stability functions
1046           rn_psi3m = -0.58
1047        CASE( 2 )                  ! Canuto A stability functions
1048           rn_psi3m = -0.64
1049        CASE( 3 )                  ! Canuto B stability functions
1050           rn_psi3m = -0.64        ! caution : constant not identified
1051        END SELECT
1052      !
1053      CASE( 3 )               ! gen
1054      !
1055      IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
1056      zpp       = 2.
1057      zmm       = 1.
1058      znn       = -0.67
1059      rn_sc_tke = 0.8
1060      rn_sc_psi = 1.07
1061      rn_psi1 = 1.
1062      rn_psi3p = 1.
1063      rn_psi2  = 1.22
1064      !
1065        SELECT CASE ( nn_stab_func )
1066        !
1067        CASE( 0, 1 )               ! G88 or KC stability functions
1068           rn_psi3m = 0.1
1069        CASE( 2 )                  ! Canuto A stability functions
1070           rn_psi3m = 0.05
1071        CASE( 3 )                  ! Canuto B stability functions
1072           rn_psi3m = 0.05         ! caution : constant not identified
1073        END SELECT
1074      !
1075      END SELECT
1076
1077      ! Initialisation of the parameters of the stability functions
1078      ! -----------------------------------------------------------
1079      !
1080      SELECT CASE ( nn_stab_func )
1081      !
1082      CASE ( 0 )             ! Galperin stability functions
1083      !
1084      IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
1085      rn_c2 = 0.
1086      rn_c3 = 0.
1087      rn_c_diff = 1.
1088      rn_c0     = 0.5544
1089      rn_cm_sf = 0.9884
1090      rn_ghmin = -0.28
1091      rn_gh0   = 0.0233
1092      rn_ghcri = 0.02
1093      !
1094      CASE ( 1 )             ! Kantha-Clayson stability functions
1095      !
1096      IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
1097      rn_c2 = 0.7
1098      rn_c3 = 0.2
1099      rn_c_diff = 1.
1100      rn_c0     = 0.5544
1101      rn_cm_sf = 0.9884
1102      rn_ghmin = -0.28
1103      rn_gh0   = 0.0233
1104      rn_ghcri = 0.02
1105      !
1106      CASE ( 2 )             ! Canuto A stability functions
1107      !
1108      IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
1109      rn_s0 = 1.5*rn_l1*rn_l5**2.
1110      rn_s1 = -rn_l4*(rn_l6+rn_l7) + 2.*rn_l4*rn_l5*(rn_l1-(1./3.)*rn_l2-rn_l3)+1.5*rn_l1*rn_l5*rn_l8
1111      rn_s2 = -(3./8.)*rn_l1*(rn_l6**2.-rn_l7**2.)
1112      rn_s4 = 2.*rn_l5
1113      rn_s5 = 2.*rn_l4
1114      rn_s6 = (2./3.)*rn_l5*(3.*rn_l3**2.-rn_l2**2.)-0.5*rn_l5*rn_l1*(3.*rn_l3-rn_l2)+0.75*rn_l1*(rn_l6-rn_l7)
1115      rn_d0 = 3*rn_l5**2.
1116      rn_d1 = rn_l5*(7.*rn_l4+3.*rn_l8)
1117      rn_d2 = rn_l5**2.*(3.*rn_l3**2.-rn_l2**2.)-0.75*(rn_l6**2.-rn_l7**2.)
1118      rn_d3 = rn_l4*(4.*rn_l4+3.*rn_l8)
1119      rn_d4 = rn_l4*(rn_l2*rn_l6-3.*rn_l3*rn_l7-rn_l5*(rn_l2**2.-rn_l3**2.))+rn_l5*rn_l8*(3.*rn_l3**2.-rn_l2**2.)
1120      rn_d5 = 0.25*(rn_l2**2.-3.*rn_l3**2.)*(rn_l6**2.-rn_l7**2.)
1121      rn_c0 = 0.5268
1122      rn_f6 = 8. / (rn_c0**6.)
1123      rn_c_diff = SQRT(2.)/(rn_c0**3.)
1124      rn_cm_sf = 0.7310
1125      rn_ghmin = -0.28
1126      rn_gh0   = 0.0329
1127      rn_ghcri = 0.03
1128      !
1129      CASE ( 3 )             ! Canuto B stability functions
1130      !
1131      IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
1132      rn_s0 = 1.5*rn_m1*rn_m5**2.
1133      rn_s1 = -rn_m4*(rn_m6+rn_m7) + 2.*rn_m4*rn_m5*(rn_m1-(1./3.)*rn_m2-rn_m3)+1.5*rn_m1*rn_m5*rn_m8
1134      rn_s2 = -(3./8.)*rn_m1*(rn_m6**2.-rn_m7**2.)
1135      rn_s4 = 2.*rn_m5
1136      rn_s5 = 2.*rn_m4
1137      rn_s6 = (2./3.)*rn_m5*(3.*rn_m3**2.-rn_m2**2.)-0.5*rn_m5*rn_m1*(3.*rn_m3-rn_m2)+0.75*rn_m1*(rn_m6-rn_m7)
1138      rn_d0 = 3*rn_m5**2.
1139      rn_d1 = rn_m5*(7.*rn_m4+3.*rn_m8)
1140      rn_d2 = rn_m5**2.*(3.*rn_m3**2.-rn_m2**2.)-0.75*(rn_m6**2.-rn_m7**2.)
1141      rn_d3 = rn_m4*(4.*rn_m4+3.*rn_m8)
1142      rn_d4 = rn_m4*(rn_m2*rn_m6-3.*rn_m3*rn_m7-rn_m5*(rn_m2**2.-rn_m3**2.))+rn_m5*rn_m8*(3.*rn_m3**2.-rn_m2**2.)
1143      rn_d5 = 0.25*(rn_m2**2.-3.*rn_m3**2.)*(rn_m6**2.-rn_m7**2.)
1144      rn_c0 = 0.5268  !!       rn_c0 = 0.5540 (Warner ...) to verify !
1145      rn_f6 = 8. / (rn_c0**6.)
1146      rn_c_diff = SQRT(2.)/(rn_c0**3.)
1147      rn_cm_sf = 0.7470
1148      rn_ghmin = -0.28
1149      rn_gh0   = 0.0444
1150      rn_ghcri = 0.0414
1151      !
1152      END SELECT
1153   
1154      ! Set Schmidt number for psi diffusion
1155      ! In the wave breaking case
1156      ! See equation 13 of Carniel et al, Ocean modelling, 30, 225-239, 2009
1157      ! or equation (17) of Burchard, JPO, 31, 3133-3145, 2001
1158      IF ((ln_sigpsi).AND.(ln_crban)) THEN
1159         zcr = SQRT(1.5*rn_sc_tke) * rn_cm_sf /vkarmn
1160         rn_sc_psi0 = vkarmn**2/(rn_psi2*rn_cm_sf**2) *              & 
1161        &             (  znn**2 - 4./3.*zcr*znn*zmm - 1./3.*zcr*znn  &
1162        &              + 2./9.*zmm*zcr**2 + 4./9.*zcr**2*zmm**2)                                 
1163      ELSE
1164         rn_sc_psi0 = rn_sc_psi
1165      ENDIF
1166 
1167      ! Shear free turbulence parameters:
1168      !
1169      rn_a_sf  = -4.*znn*SQRT(rn_sc_tke) / ( (1.+4.*zmm)*SQRT(rn_sc_tke) &
1170               &                              - SQRT(rn_sc_tke + 24.*rn_sc_psi0*rn_psi2 ) )
1171      rn_l_sf  = rn_c0 * SQRT(rn_c0/rn_cm_sf) * SQRT( ( (1. + 4.*zmm + 8.*zmm**2)*rn_sc_tke       &
1172               &                                          + 12. * rn_sc_psi0*rn_psi2 - (1. + 4.*zmm) &
1173               &                                          *SQRT(rn_sc_tke*(rn_sc_tke                &
1174               &                                                + 24.*rn_sc_psi0*rn_psi2)) )         &
1175               &                                          /(12.*znn**2.) &
1176               &                                       )
1177
1178      ! Control print
1179      !
1180      IF(lwp) THEN
1181         WRITE(numout,*)
1182         WRITE(numout,*) 'Limit values'
1183         WRITE(numout,*) '~~~~~~~~~~~~'
1184         WRITE(numout,*) 'Parameter  m = ',zmm
1185         WRITE(numout,*) 'Parameter  n = ',znn
1186         WRITE(numout,*) 'Parameter  p = ',zpp
1187         WRITE(numout,*) 'rn_psi1   = ',rn_psi1
1188         WRITE(numout,*) 'rn_psi2   = ',rn_psi2
1189         WRITE(numout,*) 'rn_psi3m  = ',rn_psi3m
1190         WRITE(numout,*) 'rn_psi3p  = ',rn_psi3p
1191         WRITE(numout,*) 'rn_sc_tke = ',rn_sc_tke
1192         WRITE(numout,*) 'rn_sc_psi = ',rn_sc_psi
1193         WRITE(numout,*) 'rn_sc_psi0 = ',rn_sc_psi0
1194         WRITE(numout,*) 'rn_c0     = ',rn_c0
1195         WRITE(numout,*)
1196         WRITE(numout,*) 'Shear free turbulence parameters:'
1197         WRITE(numout,*) 'rn_cm_sf  = ',rn_cm_sf
1198         WRITE(numout,*) 'rn_a_sf   = ',rn_a_sf
1199         WRITE(numout,*) 'rn_l_sf   = ',rn_l_sf
1200         WRITE(numout,*)
1201      ENDIF
1202
1203      ! Constants initialization
1204      rn_c02r = 1. / rn_c0**2.
1205      rn_c02  = rn_c0**2._wp
1206      rn_c03  = rn_c0**3._wp
1207      rn_c04  = rn_c0**4._wp
1208      rn_c03_sqrt2_galp = rn_c03 / SQRT(2._wp) / clim_galp
1209      rn_sbc_mb   = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg
1210      rn_sbc_std  = 3.75                           ! Surf. bound. cond. standard (prod=diss)
1211      rn_sbc_tke1 = (-rn_sc_tke*rn_crban/(rn_cm_sf*rn_a_sf*rn_l_sf))**(2./3.) ! k_eps = 53.  Dirichlet + Wave breaking
1212      rn_sbc_tke2 = 0.5 / rau0
1213      rn_sbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking
1214      rn_sbc_zs   = rn_charn/grav                                                          ! Charnock formula
1215      rn_sbc_psi1 = rn_c0**zpp * rn_sbc_tke1**zmm * rn_l_sf**znn                           ! Dirichlet + Wave breaking
1216      rn_sbc_psi2 = -0.5 * rdt * rn_c0**zpp * znn * vkarmn**znn / rn_sc_psi                   ! Neumann + NO Wave breaking
1217      rn_sbc_psi3 = -0.5 * rdt * rn_c0**zpp * rn_l_sf**znn / rn_sc_psi  * (znn + zmm*rn_a_sf) ! Neumann + Wave breaking
1218      zfact_tke = -0.5 / rn_sc_tke * rdt               ! Cst used for the Diffusion term of tke
1219      zfact_psi = -0.5 / rn_sc_psi * rdt               ! Cst used for the Diffusion term of tke
1220
1221      ! Wall proximity function
1222      zwall (:,:,:) = 1._wp * tmask(:,:,:)
1223
1224      !                               !* set vertical eddy coef. to the background value
1225      DO jk = 1, jpk
1226         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1227         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
1228         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1229         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1230      END DO
1231      !                               !* read or initialize all required files
1232      CALL gls_rst( nit000, 'READ' )
1233      !
1234   END SUBROUTINE zdf_gls_init
1235
1236   SUBROUTINE gls_rst( kt, cdrw )
1237     !!---------------------------------------------------------------------
1238     !!                   ***  ROUTINE ts_rst  ***
1239     !!                     
1240     !! ** Purpose :   Read or write TKE file (en) in restart file
1241     !!
1242     !! ** Method  :   use of IOM library
1243     !!                if the restart does not contain TKE, en is either
1244     !!                set to rn_emin or recomputed (nn_igls/=0)
1245     !!----------------------------------------------------------------------
1246     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1247     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1248     !
1249     INTEGER ::   jit, jk   ! dummy loop indices
1250     INTEGER ::   id1, id2, id3, id4, id5, id6, id7, id8
1251     INTEGER ::   ji, jj, ikbu, ikbv, ikbum1, ikbvm1
1252     REAL(wp)::   cbx, cby
1253     !!----------------------------------------------------------------------
1254     !
1255     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1256        !                                   ! ---------------
1257        IF( ln_rstart ) THEN                   !* Read the restart file
1258           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
1259           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
1260           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
1261           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
1262           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
1263           id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
1264           id7 = iom_varid( numror, 'wbotu', ldstop = .FALSE. )
1265           id8 = iom_varid( numror, 'wbotv', ldstop = .FALSE. )
1266           !
1267           IF( MIN( id1, id2, id3, id4, id5, id6, id7, id8 ) > 0 ) THEN        ! all required arrays exist
1268              CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
1269              CALL iom_get( numror, jpdom_autoglo, 'avt'   , avt    )
1270              CALL iom_get( numror, jpdom_autoglo, 'avm'   , avm    )
1271              CALL iom_get( numror, jpdom_autoglo, 'avmu'  , avmu   )
1272              CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   )
1273              CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
1274              CALL iom_get( numror, jpdom_autoglo, 'wbotu' , wbotu  )
1275              CALL iom_get( numror, jpdom_autoglo, 'wbotv' , wbotv  )
1276           ELSE                       
1277              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
1278              IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated' 
1279              en  (:,:,:) = rn_emin
1280              mxln(:,:,:) = 0.001       
1281              ! Initialize bottom stresses
1282              DO jj = 2, jpjm1
1283                DO ji = fs_2, fs_jpim1   ! vector opt.
1284                  ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
1285                  ikbum1 = MAX( ikbu-1, 1 )
1286                  ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
1287                  ikbvm1 = MAX( ikbv-1, 1 )
1288                  cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu)
1289                  cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv)
1290                  wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1)
1291                  wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1)
1292                END DO
1293              END DO
1294              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
1295           ENDIF
1296        ELSE                                   !* Start from rest
1297           IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
1298           IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated'
1299           en  (:,:,:) = rn_emin
1300           mxln(:,:,:) = 0.001       
1301           ! Initialize bottom stresses
1302           DO jj = 2, jpjm1
1303             DO ji = fs_2, fs_jpim1   ! vector opt.
1304               ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
1305               ikbum1 = MAX( ikbu-1, 1 )
1306               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
1307               ikbvm1 = MAX( ikbv-1, 1 )
1308               cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu)
1309               cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv)
1310               wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1)
1311               wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1)
1312             END DO
1313           END DO
1314        ENDIF
1315        !
1316     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1317        !                                   ! -------------------
1318        IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1319        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
1320        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt   )
1321        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm   )
1322        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu  )
1323        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  )
1324        CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln  )
1325        !
1326     ENDIF
1327     !
1328   END SUBROUTINE gls_rst
1329
1330#else
1331   !!----------------------------------------------------------------------
1332   !!   Dummy module :                                        NO TKE scheme
1333   !!----------------------------------------------------------------------
1334   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfgls = .FALSE.   !: TKE flag
1335CONTAINS
1336   SUBROUTINE zdf_gls( kt )          ! Empty routine
1337      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
1338   END SUBROUTINE zdf_gls
1339#endif
1340
1341   !!======================================================================
1342END MODULE zdfgls
Note: See TracBrowser for help on using the repository browser.