source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfgls.F90 @ 13268

Last change on this file since 13268 was 13268, checked in by clem, 3 months ago

should solve ticket #2435

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