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 NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfgls.F90 @ 13492

Last change on this file since 13492 was 13492, checked in by smasson, 4 years ago

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