source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfgls.F90 @ 12616

Last change on this file since 12616 was 12616, checked in by techene, 11 months ago

all: add e3 substitute (sometimes it requires to add ze3t/u/v/w) and limit precompiled files lines to about 130 character, OCE/ASM/asminc.F90, OCE/DOM/domzgr_substitute.h90, OCE/ISF/isfcpl.F90, OCE/SBC/sbcice_cice, OCE/CRS/crsini.F90 : add key_LF

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