source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 11442

Last change on this file since 11442 was 11442, checked in by mattmartin, 21 months ago

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

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