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 branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 8568

Last change on this file since 8568 was 8568, checked in by gm, 7 years ago

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