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

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

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/ZDF/zdfgls.F90 @ 10297

Last change on this file since 10297 was 10297, checked in by smasson, 5 years ago

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