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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 4470

Last change on this file since 4470 was 4470, checked in by trackstand2, 10 years ago

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