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

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 3798

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