New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfgls.F90 in NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/ZDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/ZDF/zdfgls.F90 @ 15717

Last change on this file since 15717 was 14922, checked in by hadcv, 3 years ago

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