source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 5021

Last change on this file since 5021 was 4839, checked in by timgraham, 6 years ago

Corrected bug fix for unitialised variables in GLS scheme when using AMM12 with reduced restarts - see ticket #1402.
Also set nn_istate=0 for AMM12 as uses too much disk space in SETTE tests.

  • Property svn:keywords set to Id
File size: 61.3 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          ! =T use Craig and Banner scheme
55   LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
56   LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
57   INTEGER  ::   nn_tkebc_surf     ! TKE surface boundary condition (=0/1)
58   INTEGER  ::   nn_tkebc_bot      ! TKE bottom boundary condition (=0/1)
59   INTEGER  ::   nn_psibc_surf     ! PSI surface boundary condition (=0/1)
60   INTEGER  ::   nn_psibc_bot      ! PSI bottom boundary condition (=0/1)
61   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
62   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
63   REAL(wp) ::   rn_clim_galp      ! Holt 2008 value for k-eps: 0.267
64   REAL(wp) ::   rn_epsmin         ! minimum value of dissipation (m2/s3)
65   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
66   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
67   REAL(wp) ::   rn_crban          ! 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      INTEGER ::   ios   ! Local integer output status for namelist read
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_ref )              ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
944      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
945901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in reference namelist', lwp )
946
947      REWIND( numnam_cfg )              ! Namelist namzdf_gls in configuration namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
948      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
949902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp )
950      IF(lwm) WRITE ( numond, namzdf_gls )
951
952      IF(lwp) THEN                     !* Control print
953         WRITE(numout,*)
954         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
955         WRITE(numout,*) '~~~~~~~~~~~~'
956         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
957         WRITE(numout,*) '      minimum value of en                           rn_emin       = ', rn_emin
958         WRITE(numout,*) '      minimum value of eps                          rn_epsmin     = ', rn_epsmin
959         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim = ', ln_length_lim
960         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp  = ', rn_clim_galp
961         WRITE(numout,*) '      TKE Surface boundary condition                nn_tkebc_surf = ', nn_tkebc_surf
962         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_tkebc_bot  = ', nn_tkebc_bot
963         WRITE(numout,*) '      PSI Surface boundary condition                nn_psibc_surf = ', nn_psibc_surf
964         WRITE(numout,*) '      PSI Bottom boundary condition                 nn_psibc_bot  = ', nn_psibc_bot
965         WRITE(numout,*) '      Craig and Banner scheme                       ln_crban      = ', ln_crban
966         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi     = ', ln_sigpsi
967         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
968         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
969         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
970         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
971         WRITE(numout,*) '   Hard coded parameters'
972         WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro
973         WRITE(numout,*) '      Bottom roughness (m)                          hbro          = ', hbro
974      ENDIF
975
976      !                                !* allocate gls arrays
977      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
978
979      !                                !* Check of some namelist values
980      IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' )
981      IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' )
982      IF( nn_tkebc_bot  < 0 .OR. nn_tkebc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' )
983      IF( nn_psibc_bot  < 0 .OR. nn_psibc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' )
984      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
985      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
986
987      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
988      !
989      CASE( 0 )                              ! k-kl  (Mellor-Yamada)
990         !
991         IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
992         rpp     = 0._wp
993         rmm     = 1._wp
994         rnn     = 1._wp
995         rsc_tke = 1.96_wp
996         rsc_psi = 1.96_wp
997         rpsi1   = 0.9_wp
998         rpsi3p  = 1._wp
999         rpsi2   = 0.5_wp
1000         !
1001         SELECT CASE ( nn_stab_func )
1002         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
1003         CASE( 2 )      ;   rpsi3m = 2.38_wp       ! Canuto A stability functions
1004         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
1005         END SELECT
1006         !
1007      CASE( 1 )                              ! k-eps
1008         !
1009         IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
1010         rpp     =  3._wp
1011         rmm     =  1.5_wp
1012         rnn     = -1._wp
1013         rsc_tke =  1._wp
1014         rsc_psi =  1.3_wp  ! Schmidt number for psi
1015         rpsi1   =  1.44_wp
1016         rpsi3p  =  1._wp
1017         rpsi2   =  1.92_wp
1018         !
1019         SELECT CASE ( nn_stab_func )
1020         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
1021         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
1022         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
1023         END SELECT
1024         !
1025      CASE( 2 )                              ! k-omega
1026         !
1027         IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
1028         rpp     = -1._wp
1029         rmm     =  0.5_wp
1030         rnn     = -1._wp
1031         rsc_tke =  2._wp
1032         rsc_psi =  2._wp
1033         rpsi1   =  0.555_wp
1034         rpsi3p  =  1._wp
1035         rpsi2   =  0.833_wp
1036         !
1037         SELECT CASE ( nn_stab_func )
1038         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
1039         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
1040         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
1041         END SELECT
1042         !
1043      CASE( 3 )                              ! generic
1044         !
1045         IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
1046         rpp     = 2._wp
1047         rmm     = 1._wp
1048         rnn     = -0.67_wp
1049         rsc_tke = 0.8_wp
1050         rsc_psi = 1.07_wp
1051         rpsi1   = 1._wp
1052         rpsi3p  = 1._wp
1053         rpsi2   = 1.22_wp
1054         !
1055         SELECT CASE ( nn_stab_func )
1056         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
1057         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
1058         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
1059         END SELECT
1060         !
1061      END SELECT
1062
1063      !
1064      SELECT CASE ( nn_stab_func )     !* set the parameters of the stability functions
1065      !
1066      CASE ( 0 )                             ! Galperin stability functions
1067         !
1068         IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
1069         rc2     =  0._wp
1070         rc3     =  0._wp
1071         rc_diff =  1._wp
1072         rc0     =  0.5544_wp
1073         rcm_sf  =  0.9884_wp
1074         rghmin  = -0.28_wp
1075         rgh0    =  0.0233_wp
1076         rghcri  =  0.02_wp
1077         !
1078      CASE ( 1 )                             ! Kantha-Clayson stability functions
1079         !
1080         IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
1081         rc2     =  0.7_wp
1082         rc3     =  0.2_wp
1083         rc_diff =  1._wp
1084         rc0     =  0.5544_wp
1085         rcm_sf  =  0.9884_wp
1086         rghmin  = -0.28_wp
1087         rgh0    =  0.0233_wp
1088         rghcri  =  0.02_wp
1089         !
1090      CASE ( 2 )                             ! Canuto A stability functions
1091         !
1092         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
1093         rs0 = 1.5_wp * rl1 * rl5*rl5
1094         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
1095         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
1096         rs4 = 2._wp * rl5
1097         rs5 = 2._wp * rl4
1098         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
1099            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
1100         rd0 = 3._wp * rl5*rl5
1101         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
1102         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
1103         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
1104         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
1105         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
1106         rc0 = 0.5268_wp
1107         rf6 = 8._wp / (rc0**6._wp)
1108         rc_diff = SQRT(2._wp) / (rc0**3._wp)
1109         rcm_sf  =  0.7310_wp
1110         rghmin  = -0.28_wp
1111         rgh0    =  0.0329_wp
1112         rghcri  =  0.03_wp
1113         !
1114      CASE ( 3 )                             ! Canuto B stability functions
1115         !
1116         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
1117         rs0 = 1.5_wp * rm1 * rm5*rm5
1118         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
1119         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
1120         rs4 = 2._wp * rm5
1121         rs5 = 2._wp * rm4
1122         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)
1123         rd0 = 3._wp * rm5*rm5
1124         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
1125         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
1126         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
1127         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
1128         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
1129         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
1130         rf6 = 8._wp / ( rc0**6._wp )
1131         rc_diff = SQRT(2._wp)/(rc0**3.)
1132         rcm_sf  =  0.7470_wp
1133         rghmin  = -0.28_wp
1134         rgh0    =  0.0444_wp
1135         rghcri  =  0.0414_wp
1136         !
1137      END SELECT
1138   
1139      !                                !* Set Schmidt number for psi diffusion in the wave breaking case
1140      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1141      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1142      IF( ln_sigpsi .AND. ln_crban ) THEN
1143         zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn
1144         rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf )                       & 
1145        &         * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn   &
1146        &           + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm )                                 
1147      ELSE
1148         rsc_psi0 = rsc_psi
1149      ENDIF
1150 
1151      !                                !* Shear free turbulence parameters
1152      !
1153      ra_sf  = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke )   &
1154         &                                      - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
1155      rl_sf  = rc0 * SQRT( rc0 / rcm_sf )                                                                   &
1156         &         * SQRT(  (  (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke                                &
1157         &                   + 12._wp * rsc_psi0 * rpsi2                                                    &
1158         &                   - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) )  )    &
1159         &                / ( 12._wp*rnn*rnn )                                                              )
1160
1161      !
1162      IF(lwp) THEN                     !* Control print
1163         WRITE(numout,*)
1164         WRITE(numout,*) 'Limit values'
1165         WRITE(numout,*) '~~~~~~~~~~~~'
1166         WRITE(numout,*) 'Parameter  m = ',rmm
1167         WRITE(numout,*) 'Parameter  n = ',rnn
1168         WRITE(numout,*) 'Parameter  p = ',rpp
1169         WRITE(numout,*) 'rpsi1   = ',rpsi1
1170         WRITE(numout,*) 'rpsi2   = ',rpsi2
1171         WRITE(numout,*) 'rpsi3m  = ',rpsi3m
1172         WRITE(numout,*) 'rpsi3p  = ',rpsi3p
1173         WRITE(numout,*) 'rsc_tke = ',rsc_tke
1174         WRITE(numout,*) 'rsc_psi = ',rsc_psi
1175         WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0
1176         WRITE(numout,*) 'rc0     = ',rc0
1177         WRITE(numout,*)
1178         WRITE(numout,*) 'Shear free turbulence parameters:'
1179         WRITE(numout,*) 'rcm_sf  = ',rcm_sf
1180         WRITE(numout,*) 'ra_sf   = ',ra_sf
1181         WRITE(numout,*) 'rl_sf   = ',rl_sf
1182         WRITE(numout,*)
1183      ENDIF
1184
1185      !                                !* Constants initialization
1186      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
1187      rc03  = rc02 * rc0
1188      rc04  = rc03 * rc0
1189      rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp
1190      rsbc_mb   = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp)               ! Surf. bound. cond. from Mellor and Blumberg
1191      rsbc_std  = 3.75_wp                                                  ! Surf. bound. cond. standard (prod=diss)
1192      rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.  Dirichlet + Wave breaking
1193      rsbc_tke2 = 0.5_wp / rau0
1194      rsbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking
1195      rsbc_zs   = rn_charn / grav                                                        ! Charnock formula
1196      rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn                           ! Dirichlet + Wave breaking
1197      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking
1198      rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + Wave breaking
1199      rfact_tke = -0.5_wp / rsc_tke * rdt               ! Cst used for the Diffusion term of tke
1200      rfact_psi = -0.5_wp / rsc_psi * rdt               ! Cst used for the Diffusion term of tke
1201
1202      !                                !* Wall proximity function
1203      zwall (:,:,:) = 1._wp * tmask(:,:,:)
1204
1205      !                                !* set vertical eddy coef. to the background value
1206      DO jk = 1, jpk
1207         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1208         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
1209         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1210         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1211      END DO
1212      !                             
1213      CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
1214      !
1215      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init')
1216      !
1217   END SUBROUTINE zdf_gls_init
1218
1219
1220   SUBROUTINE gls_rst( kt, cdrw )
1221      !!---------------------------------------------------------------------
1222      !!                   ***  ROUTINE ts_rst  ***
1223      !!                     
1224      !! ** Purpose :   Read or write TKE file (en) in restart file
1225      !!
1226      !! ** Method  :   use of IOM library
1227      !!                if the restart does not contain TKE, en is either
1228      !!                set to rn_emin or recomputed (nn_igls/=0)
1229      !!----------------------------------------------------------------------
1230      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1231      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1232      !
1233      INTEGER ::   jit, jk   ! dummy loop indices
1234      INTEGER ::   id1, id2, id3, id4, id5, id6
1235      INTEGER ::   ji, jj, ikbu, ikbv
1236      REAL(wp)::   cbx, cby
1237      !!----------------------------------------------------------------------
1238      !
1239      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1240         !                                   ! ---------------
1241         IF( ln_rstart ) THEN                   !* Read the restart file
1242            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
1243            id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
1244            id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
1245            id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
1246            id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
1247            id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
1248            !
1249            IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
1250               CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
1251               CALL iom_get( numror, jpdom_autoglo, 'avt'   , avt    )
1252               CALL iom_get( numror, jpdom_autoglo, 'avm'   , avm    )
1253               CALL iom_get( numror, jpdom_autoglo, 'avmu'  , avmu   )
1254               CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   )
1255               CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
1256            ELSE                       
1257               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
1258               en  (:,:,:) = rn_emin
1259               mxln(:,:,:) = 0.001       
1260               avt_k (:,:,:) = avt (:,:,:)
1261               avm_k (:,:,:) = avm (:,:,:)
1262               avmu_k(:,:,:) = avmu(:,:,:)
1263               avmv_k(:,:,:) = avmv(:,:,:)
1264               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
1265            ENDIF
1266         ELSE                                   !* Start from rest
1267            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
1268            en  (:,:,:) = rn_emin
1269            mxln(:,:,:) = 0.001       
1270         ENDIF
1271         !
1272      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1273         !                                   ! -------------------
1274         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1275         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
1276         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
1277         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
1278         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
1279         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
1280         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   )
1281         !
1282      ENDIF
1283      !
1284   END SUBROUTINE gls_rst
1285
1286#else
1287   !!----------------------------------------------------------------------
1288   !!   Dummy module :                                        NO TKE scheme
1289   !!----------------------------------------------------------------------
1290   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfgls = .FALSE.   !: TKE flag
1291CONTAINS
1292   SUBROUTINE zdf_gls_init           ! Empty routine
1293      WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
1294   END SUBROUTINE zdf_gls_init
1295   SUBROUTINE zdf_gls( kt )          ! Empty routine
1296      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
1297   END SUBROUTINE zdf_gls
1298   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
1299      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1300      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1301      WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw
1302   END SUBROUTINE gls_rst
1303#endif
1304
1305   !!======================================================================
1306END MODULE zdfgls
1307
Note: See TracBrowser for help on using the repository browser.