source: NEMO/branches/2020/dev_12905_xios_restart/src/OCE/ZDF/zdfgls.F90 @ 12950

Last change on this file since 12950 was 12950, checked in by andmirek, 6 months ago

Ticket #2462: new XIOS restart read/write interfaces

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