source: NEMO/trunk/src/OCE/ZDF/zdfgls.F90 @ 13237

Last change on this file since 13237 was 13237, checked in by smasson, 3 months ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

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