New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfgls.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

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