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

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 3804

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

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