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/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdfgls.F90 @ 13249

Last change on this file since 13249 was 13249, checked in by clem, 4 years ago

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