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

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

Last change on this file since 4677 was 4624, checked in by acc, 10 years ago

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