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/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/ZDF/zdfgls.F90 @ 12210

Last change on this file since 12210 was 12210, checked in by cetlod, 4 years ago

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