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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 2412

Last change on this file since 2412 was 2409, checked in by cetlod, 14 years ago

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