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

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 7990

Last change on this file since 7990 was 7990, checked in by gm, 7 years ago

#1880 (HPC-09): OPA remove avmu, avmv from zdf modules + move CALL tke(gls)_rst & gls_rst in zdftke(gls) + rename zdftmx and zdfqiao

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