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

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

source: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/ZDF/zdfgls.F90 @ 10455

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

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 59.2 KB
Line 
1MODULE zdfgls
2   !!======================================================================
3   !!                       ***  MODULE  zdfgls  ***
4   !! Ocean physics:  vertical mixing coefficient computed from the gls
5   !!                 turbulent closure parameterization
6   !!======================================================================
7   !! History :  3.0  !  2009-09  (G. Reffray)  Original code
8   !!            3.3  !  2010-10  (C. Bricaud)  Add in the reference
9   !!            4.0  !  2017-04  (G. Madec)  remove CPP keys & avm at t-point only
10   !!             -   !  2017-05  (G. Madec)  add top friction as boundary condition
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   zdf_gls       : update momentum and tracer Kz from a gls scheme
15   !!   zdf_gls_init  : initialization, namelist read, and parameters control
16   !!   gls_rst       : read/write gls restart in ocean restart file
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and active tracers
19   USE dom_oce        ! ocean space and time domain
20   USE domvvl         ! ocean space and time domain : variable volume layer
21   USE zdfdrg  , ONLY : 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      CALL mpp_sum ( 'zdfgls', zdf_gls_alloc )
122      IF( zdf_gls_alloc /= 0 )   CALL ctl_stop( 'STOP', '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, in fact not used for jk = 2 (see surface conditions)
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, in fact not used for jk = ibotm1 (see bottom conditions)
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      ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
795      zstm(:,:,jpk) = 0. 
796      DO jj = 2, jpjm1                ! update bottom with good values
797         DO ji = fs_2, fs_jpim1   ! vector opt.
798            zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
799         END DO
800      END DO
801
802      zstt(:,:,  1) = wmask(:,:,  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
803      zstt(:,:,jpk) = wmask(:,:,jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
804
805!!gm should be done for ISF (top boundary cond.)
806!!gm so, totally new staff needed!!gm
807
808      ! Compute diffusivities/viscosities
809      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
810      !  -> yes BUT p_avm(:,:1) and p_avm(:,:jpk) are used when we compute zd_lw(:,:2) and zd_up(:,:jpkm1). These values are
811      !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
812      !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
813      DO jk = 1, jpk
814         DO jj = 2, jpjm1
815            DO ji = fs_2, fs_jpim1   ! vector opt.
816               zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
817               zavt  = zsqen * zstt(ji,jj,jk)
818               zavm  = zsqen * zstm(ji,jj,jk)
819               p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
820               p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
821            END DO
822         END DO
823      END DO
824      p_avt(:,:,1) = 0._wp
825      !
826      IF(ln_ctl) THEN
827         CALL prt_ctl( tab3d_1=en   , clinfo1=' gls  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
828         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls  - m: ', kdim=jpk )
829      ENDIF
830      !
831   END SUBROUTINE zdf_gls
832
833
834   SUBROUTINE zdf_gls_init
835      !!----------------------------------------------------------------------
836      !!                  ***  ROUTINE zdf_gls_init  ***
837      !!                     
838      !! ** Purpose :   Initialization of the vertical eddy diffivity and
839      !!              viscosity computed using a GLS turbulent closure scheme
840      !!
841      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
842      !!
843      !! ** input   :   Namlist namzdf_gls
844      !!
845      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
846      !!
847      !!----------------------------------------------------------------------
848      INTEGER ::   jk    ! dummy loop indices
849      INTEGER ::   ios   ! Local integer output status for namelist read
850      REAL(wp)::   zcr   ! local scalar
851      !!
852      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
853         &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
854         &            rn_crban, rn_charn, rn_frac_hs,        &
855         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
856         &            nn_stab_func, nn_clos
857      !!----------------------------------------------------------
858      !
859      REWIND( numnam_ref )              ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
860      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
861901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_gls in reference namelist', lwp )
862
863      REWIND( numnam_cfg )              ! Namelist namzdf_gls in configuration namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
864      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
865902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp )
866      IF(lwm) WRITE ( numond, namzdf_gls )
867
868      IF(lwp) THEN                     !* Control print
869         WRITE(numout,*)
870         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
871         WRITE(numout,*) '~~~~~~~~~~~~'
872         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
873         WRITE(numout,*) '      minimum value of en                           rn_emin        = ', rn_emin
874         WRITE(numout,*) '      minimum value of eps                          rn_epsmin      = ', rn_epsmin
875         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim  = ', ln_length_lim
876         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp   = ', rn_clim_galp
877         WRITE(numout,*) '      TKE Surface boundary condition                nn_bc_surf     = ', nn_bc_surf
878         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
879         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
880         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
881         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
882         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
883         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs
884         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
885         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
886         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
887         WRITE(numout,*)
888         WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
889         WRITE(numout,*) '      top    ocean cavity roughness (m)             rn_z0(_top)   = ', r_z0_top
890         WRITE(numout,*) '      Bottom seafloor     roughness (m)             rn_z0(_bot)   = ', r_z0_bot
891         WRITE(numout,*)
892      ENDIF
893
894      !                                !* allocate GLS arrays
895      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
896
897      !                                !* Check of some namelist values
898      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' )
899      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' )
900      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' )
901      IF( nn_z0_met == 3  .AND. .NOT.ln_wave     )   CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )
902      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' )
903      IF( nn_clos      < 0 .OR. nn_clos      > 3 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
904
905      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
906      !
907      CASE( 0 )                              ! k-kl  (Mellor-Yamada)
908         !
909         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
910         IF(lwp) WRITE(numout,*)
911         rpp     = 0._wp
912         rmm     = 1._wp
913         rnn     = 1._wp
914         rsc_tke = 1.96_wp
915         rsc_psi = 1.96_wp
916         rpsi1   = 0.9_wp
917         rpsi3p  = 1._wp
918         rpsi2   = 0.5_wp
919         !
920         SELECT CASE ( nn_stab_func )
921         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
922         CASE( 2 )      ;   rpsi3m = 2.62_wp       ! Canuto A stability functions
923         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
924         END SELECT
925         !
926      CASE( 1 )                              ! k-eps
927         !
928         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen'
929         IF(lwp) WRITE(numout,*)
930         rpp     =  3._wp
931         rmm     =  1.5_wp
932         rnn     = -1._wp
933         rsc_tke =  1._wp
934         rsc_psi =  1.2_wp  ! Schmidt number for psi
935         rpsi1   =  1.44_wp
936         rpsi3p  =  1._wp
937         rpsi2   =  1.92_wp
938         !
939         SELECT CASE ( nn_stab_func )
940         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
941         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
942         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
943         END SELECT
944         !
945      CASE( 2 )                              ! k-omega
946         !
947         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen'
948         IF(lwp) WRITE(numout,*)
949         rpp     = -1._wp
950         rmm     =  0.5_wp
951         rnn     = -1._wp
952         rsc_tke =  2._wp
953         rsc_psi =  2._wp
954         rpsi1   =  0.555_wp
955         rpsi3p  =  1._wp
956         rpsi2   =  0.833_wp
957         !
958         SELECT CASE ( nn_stab_func )
959         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
960         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
961         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
962         END SELECT
963         !
964      CASE( 3 )                              ! generic
965         !
966         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen'
967         IF(lwp) WRITE(numout,*)
968         rpp     = 2._wp
969         rmm     = 1._wp
970         rnn     = -0.67_wp
971         rsc_tke = 0.8_wp
972         rsc_psi = 1.07_wp
973         rpsi1   = 1._wp
974         rpsi3p  = 1._wp
975         rpsi2   = 1.22_wp
976         !
977         SELECT CASE ( nn_stab_func )
978         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
979         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
980         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
981         END SELECT
982         !
983      END SELECT
984
985      !
986      SELECT CASE ( nn_stab_func )     !* set the parameters of the stability functions
987      !
988      CASE ( 0 )                             ! Galperin stability functions
989         !
990         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin'
991         rc2     =  0._wp
992         rc3     =  0._wp
993         rc_diff =  1._wp
994         rc0     =  0.5544_wp
995         rcm_sf  =  0.9884_wp
996         rghmin  = -0.28_wp
997         rgh0    =  0.0233_wp
998         rghcri  =  0.02_wp
999         !
1000      CASE ( 1 )                             ! Kantha-Clayson stability functions
1001         !
1002         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson'
1003         rc2     =  0.7_wp
1004         rc3     =  0.2_wp
1005         rc_diff =  1._wp
1006         rc0     =  0.5544_wp
1007         rcm_sf  =  0.9884_wp
1008         rghmin  = -0.28_wp
1009         rgh0    =  0.0233_wp
1010         rghcri  =  0.02_wp
1011         !
1012      CASE ( 2 )                             ! Canuto A stability functions
1013         !
1014         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A'
1015         rs0 = 1.5_wp * rl1 * rl5*rl5
1016         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
1017         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
1018         rs4 = 2._wp * rl5
1019         rs5 = 2._wp * rl4
1020         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
1021            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
1022         rd0 = 3._wp * rl5*rl5
1023         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
1024         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
1025         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
1026         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
1027         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
1028         rc0 = 0.5268_wp
1029         rf6 = 8._wp / (rc0**6._wp)
1030         rc_diff = SQRT(2._wp) / (rc0**3._wp)
1031         rcm_sf  =  0.7310_wp
1032         rghmin  = -0.28_wp
1033         rgh0    =  0.0329_wp
1034         rghcri  =  0.03_wp
1035         !
1036      CASE ( 3 )                             ! Canuto B stability functions
1037         !
1038         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B'
1039         rs0 = 1.5_wp * rm1 * rm5*rm5
1040         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
1041         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
1042         rs4 = 2._wp * rm5
1043         rs5 = 2._wp * rm4
1044         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)
1045         rd0 = 3._wp * rm5*rm5
1046         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
1047         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
1048         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
1049         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
1050         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
1051         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
1052         rf6 = 8._wp / ( rc0**6._wp )
1053         rc_diff = SQRT(2._wp)/(rc0**3.)
1054         rcm_sf  =  0.7470_wp
1055         rghmin  = -0.28_wp
1056         rgh0    =  0.0444_wp
1057         rghcri  =  0.0414_wp
1058         !
1059      END SELECT
1060   
1061      !                                !* Set Schmidt number for psi diffusion in the wave breaking case
1062      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1063      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1064      IF( ln_sigpsi ) THEN
1065         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
1066         ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
1067         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
1068         ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn
1069         rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.))
1070      ELSE
1071         rsc_psi0 = rsc_psi
1072      ENDIF
1073 
1074      !                                !* Shear free turbulence parameters
1075      !
1076      ra_sf  = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) &
1077               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
1078
1079      IF ( rn_crban==0._wp ) THEN
1080         rl_sf = vkarmn
1081      ELSE
1082         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp) * rsc_tke        &
1083            &                                            + 12._wp*rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm)   &
1084            &                                                     *SQRT(rsc_tke*(rsc_tke                 &
1085            &                                                        + 24._wp*rsc_psi0*rpsi2)) )         &
1086            &                                              /(12._wp*rnn**2.)                             )
1087      ENDIF
1088
1089      !
1090      IF(lwp) THEN                     !* Control print
1091         WRITE(numout,*)
1092         WRITE(numout,*) '   Limit values :'
1093         WRITE(numout,*) '      Parameter  m = ', rmm
1094         WRITE(numout,*) '      Parameter  n = ', rnn
1095         WRITE(numout,*) '      Parameter  p = ', rpp
1096         WRITE(numout,*) '      rpsi1    = ', rpsi1
1097         WRITE(numout,*) '      rpsi2    = ', rpsi2
1098         WRITE(numout,*) '      rpsi3m   = ', rpsi3m
1099         WRITE(numout,*) '      rpsi3p   = ', rpsi3p
1100         WRITE(numout,*) '      rsc_tke  = ', rsc_tke
1101         WRITE(numout,*) '      rsc_psi  = ', rsc_psi
1102         WRITE(numout,*) '      rsc_psi0 = ', rsc_psi0
1103         WRITE(numout,*) '      rc0      = ', rc0
1104         WRITE(numout,*)
1105         WRITE(numout,*) '   Shear free turbulence parameters:'
1106         WRITE(numout,*) '      rcm_sf   = ', rcm_sf
1107         WRITE(numout,*) '      ra_sf    = ', ra_sf
1108         WRITE(numout,*) '      rl_sf    = ', rl_sf
1109      ENDIF
1110
1111      !                                !* Constants initialization
1112      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
1113      rc03  = rc02 * rc0
1114      rc04  = rc03 * rc0
1115      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
1116      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking
1117      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
1118      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer
1119      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
1120      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
1121      rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
1122      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
1123      !
1124      rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
1125      rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke
1126      !
1127      !                                !* Wall proximity function
1128!!gm tmask or wmask ????
1129      zwall(:,:,:) = 1._wp * tmask(:,:,:)
1130
1131      !                                !* read or initialize all required files 
1132      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxl_n)
1133      !
1134      IF( lwxios ) THEN
1135         CALL iom_set_rstw_var_active('en')
1136         CALL iom_set_rstw_var_active('avt_k')
1137         CALL iom_set_rstw_var_active('avm_k')
1138         CALL iom_set_rstw_var_active('hmxl_n')
1139      ENDIF
1140      !
1141   END SUBROUTINE zdf_gls_init
1142
1143
1144   SUBROUTINE gls_rst( kt, cdrw )
1145      !!---------------------------------------------------------------------
1146      !!                   ***  ROUTINE gls_rst  ***
1147      !!                     
1148      !! ** Purpose :   Read or write TKE file (en) in restart file
1149      !!
1150      !! ** Method  :   use of IOM library
1151      !!                if the restart does not contain TKE, en is either
1152      !!                set to rn_emin or recomputed (nn_igls/=0)
1153      !!----------------------------------------------------------------------
1154      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
1155      !!
1156      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1157      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1158      !
1159      INTEGER ::   jit, jk   ! dummy loop indices
1160      INTEGER ::   id1, id2, id3, id4
1161      INTEGER ::   ji, jj, ikbu, ikbv
1162      REAL(wp)::   cbx, cby
1163      !!----------------------------------------------------------------------
1164      !
1165      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1166         !                                   ! ---------------
1167         IF( ln_rstart ) THEN                   !* Read the restart file
1168            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. )
1169            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
1170            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
1171            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
1172            !
1173            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
1174               CALL iom_get( numror, jpdom_autoglo, 'en'    , en    , ldxios = lrxios )
1175               CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k , ldxios = lrxios )
1176               CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k , ldxios = lrxios )
1177               CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n, ldxios = lrxios )
1178            ELSE                       
1179               IF(lwp) WRITE(numout,*)
1180               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values'
1181               en    (:,:,:) = rn_emin
1182               hmxl_n(:,:,:) = 0.05_wp
1183               ! avt_k, avm_k already set to the background value in zdf_phy_init
1184            ENDIF
1185         ELSE                                   !* Start from rest
1186            IF(lwp) WRITE(numout,*)
1187            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values'
1188            en    (:,:,:) = rn_emin
1189            hmxl_n(:,:,:) = 0.05_wp
1190            ! avt_k, avm_k already set to the background value in zdf_phy_init
1191         ENDIF
1192         !
1193      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1194         !                                   ! -------------------
1195         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1196         IF( lwxios ) CALL iom_swap(      cwxios_context         )
1197         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en    , ldxios = lwxios )
1198         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k , ldxios = lwxios )
1199         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k , ldxios = lwxios )
1200         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n, ldxios = lwxios )
1201         IF( lwxios ) CALL iom_swap(      cxios_context          )
1202         !
1203      ENDIF
1204      !
1205   END SUBROUTINE gls_rst
1206
1207   !!======================================================================
1208END MODULE zdfgls
1209
Note: See TracBrowser for help on using the repository browser.