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

Last change on this file was 14156, checked in by jchanut, 4 months ago

#2560, correct bottom Neumann boundary conditions and bottom friction velocities

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