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 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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