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

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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 9023

Last change on this file since 9023 was 9023, checked in by timgraham, 6 years ago

Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC
  • Property svn:keywords set to Id
File size: 58.1 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 timing         ! Timing
34   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   zdf_gls        ! called in zdfphy
40   PUBLIC   zdf_gls_init   ! called in zdfphy
41   PUBLIC   gls_rst        ! called in zdfphy
42
43   !
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxl_n    !: now mixing length
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_surf !: Squared surface velocity scale at T-points
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_top  !: Squared top     velocity scale at T-points
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_bot  !: Squared bottom  velocity scale at T-points
49
50   !                              !! ** Namelist  namzdf_gls  **
51   LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
52   LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
53   INTEGER  ::   nn_bc_surf        ! surface boundary condition (=0/1)
54   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1)
55   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation
56   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
57   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
58   REAL(wp) ::   rn_clim_galp      ! Holt 2008 value for k-eps: 0.267
59   REAL(wp) ::   rn_epsmin         ! minimum value of dissipation (m2/s3)
60   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
61   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
62   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
63   REAL(wp) ::   rn_hsro           ! Minimum surface roughness
64   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)
65
66   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
67   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
68   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn   
69   REAL(wp) ::   rghmin        = -0.28_wp
70   REAL(wp) ::   rgh0          =  0.0329_wp
71   REAL(wp) ::   rghcri        =  0.03_wp
72   REAL(wp) ::   ra1           =  0.92_wp
73   REAL(wp) ::   ra2           =  0.74_wp
74   REAL(wp) ::   rb1           = 16.60_wp
75   REAL(wp) ::   rb2           = 10.10_wp         
76   REAL(wp) ::   re2           =  1.33_wp         
77   REAL(wp) ::   rl1           =  0.107_wp
78   REAL(wp) ::   rl2           =  0.0032_wp
79   REAL(wp) ::   rl3           =  0.0864_wp
80   REAL(wp) ::   rl4           =  0.12_wp
81   REAL(wp) ::   rl5           = 11.9_wp
82   REAL(wp) ::   rl6           =  0.4_wp
83   REAL(wp) ::   rl7           =  0.0_wp
84   REAL(wp) ::   rl8           =  0.48_wp
85   REAL(wp) ::   rm1           =  0.127_wp
86   REAL(wp) ::   rm2           =  0.00336_wp
87   REAL(wp) ::   rm3           =  0.0906_wp
88   REAL(wp) ::   rm4           =  0.101_wp
89   REAL(wp) ::   rm5           = 11.2_wp
90   REAL(wp) ::   rm6           =  0.4_wp
91   REAL(wp) ::   rm7           =  0.0_wp
92   REAL(wp) ::   rm8           =  0.318_wp
93   REAL(wp) ::   rtrans        =  0.1_wp
94   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
95   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        -
96   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        -
97   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
98   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
99   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        -
100   REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5                     !     -           -           -        -
101   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        -
102   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        -
103   !
104   REAL(wp) ::   r2_3 = 2._wp/3._wp   ! constant=2/3
105
106   !! * Substitutions
107#  include "vectopt_loop_substitute.h90"
108   !!----------------------------------------------------------------------
109   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
110   !! $Id$
111   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
112   !!----------------------------------------------------------------------
113CONTAINS
114
115   INTEGER FUNCTION zdf_gls_alloc()
116      !!----------------------------------------------------------------------
117      !!                ***  FUNCTION zdf_gls_alloc  ***
118      !!----------------------------------------------------------------------
119      ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) ,                     &
120         &      zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc )
121         !
122      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
123      IF( zdf_gls_alloc /= 0 )   CALL ctl_warn('zdf_gls_alloc: failed to allocate arrays')
124   END FUNCTION zdf_gls_alloc
125
126
127   SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt )
128      !!----------------------------------------------------------------------
129      !!                   ***  ROUTINE zdf_gls  ***
130      !!
131      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
132      !!              coefficients using the GLS turbulent closure scheme.
133      !!----------------------------------------------------------------------
134      USE zdf_oce , ONLY : en, avtb, avmb   ! ocean vertical physics
135      !!
136      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step
137      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
138      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
139      !
140      INTEGER  ::   ji, jj, jk    ! dummy loop arguments
141      INTEGER  ::   ibot, ibotm1  ! local integers
142      INTEGER  ::   itop, itopp1  !   -       -
143      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1 , zex2  ! local scalars
144      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof, zdir  !   -      -
145      REAL(wp) ::   zratio, zrn2, zflxb, sh     , z_en  !   -      -
146      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
147      REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      -
148      REAL(wp) ::   zmsku, zmskv                        !   -      -
149      REAL(wp), DIMENSION(jpi,jpj)     ::   zdep
150      REAL(wp), DIMENSION(jpi,jpj)     ::   zkar
151      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves
152      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves)
153      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before
154      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before
155      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps         ! dissipation rate
156      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
157      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   psi         ! psi at time now
158      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
159      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zstt, zstm  ! stability function on tracer and momentum
160      !!--------------------------------------------------------------------
161      !
162      IF( ln_timing )   CALL timing_start('zdf_gls')
163      !
164      ! Preliminary computing
165
166      ustar2_surf(:,:) = 0._wp   ;         psi(:,:,:) = 0._wp   
167      ustar2_top (:,:) = 0._wp   ;   zwall_psi(:,:,:) = 0._wp
168      ustar2_bot (:,:) = 0._wp
169
170      ! Compute surface, top and bottom friction at T-points
171      DO jj = 2, jpjm1         
172         DO ji = fs_2, fs_jpim1   ! vector opt.         
173            !
174            ! surface friction
175            ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
176            !   
177!!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that...
178          ! bottom friction (explicit before friction)
179          zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
180          zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
181          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
182             &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
183         END DO
184      END DO
185      IF( ln_isfcav ) THEN       !top friction
186         DO jj = 2, jpjm1
187            DO ji = fs_2, fs_jpim1   ! vector opt.
188               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
189               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
190               ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
191                  &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
192            END DO
193         END DO
194      ENDIF
195   
196      SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==!
197      CASE ( 0 )                          ! Constant roughness         
198         zhsro(:,:) = rn_hsro
199      CASE ( 1 )             ! Standard Charnock formula
200         zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro )
201      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
202!!gm faster coding : the 2 comment lines should be used
203!!gm         zcof = 2._wp * 0.6_wp / 28._wp
204!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) )  )       ! Wave age (eq. 10)
205         zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) )         ! Wave age (eq. 10)
206         zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro)          ! zhsro = rn_frac_hs * Hsw (eq. 11)
207      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
208         zhsro(:,:) = rn_frac_hs * hsw(:,:)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )
209      END SELECT
210      !
211      DO jk = 2, jpkm1              !==  Compute dissipation rate  ==!
212         DO jj = 1, jpjm1
213            DO ji = 1, jpim1
214               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
215            END DO
216         END DO
217      END DO
218
219      ! Save tke at before time step
220      eb    (:,:,:) = en    (:,:,:)
221      hmxl_b(:,:,:) = hmxl_n(:,:,:)
222
223      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
224         DO jk = 2, jpkm1
225            DO jj = 2, jpjm1 
226               DO ji = fs_2, fs_jpim1   ! vector opt.
227                  zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
228                  zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
229                  zcoef = ( zup / MAX( zdown, rsmall ) )
230                  zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk)
231               END DO
232            END DO
233         END DO
234      ENDIF
235
236      !!---------------------------------!!
237      !!   Equation to prognostic k      !!
238      !!---------------------------------!!
239      !
240      ! Now Turbulent kinetic energy (output in en)
241      ! -------------------------------
242      ! Resolution of a tridiagonal linear system by a "methode de chasse"
243      ! computation from level 2 to jpkm1  (e(1) computed after and e(jpk)=0 ).
244      ! The surface boundary condition are set after
245      ! The bottom boundary condition are also set after. In standard e(bottom)=0.
246      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
247      ! Warning : after this step, en : right hand side of the matrix
248
249      DO jk = 2, jpkm1
250         DO jj = 2, jpjm1
251            DO ji = 2, jpim1
252               !
253               buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
254               !
255               diss = eps(ji,jj,jk)                         ! dissipation
256               !
257               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)
258               !
259               zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk)          ! production term
260               zdiss = zdir*(diss/en(ji,jj,jk))   +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
261!!gm better coding, identical results
262!               zesh2 =   p_sh2(ji,jj,jk) + zdir*buoy               ! production term
263!               zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term
264!!gm
265               !
266               ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
267               ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style)
268               ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
269               ! Otherwise, this should be rsc_psi/rsc_psi0
270               IF( ln_sigpsi ) THEN
271                  zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
272                  zwall_psi(ji,jj,jk) = rsc_psi /   & 
273                     &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
274               ELSE
275                  zwall_psi(ji,jj,jk) = 1._wp
276               ENDIF
277               !
278               ! building the matrix
279               zcof = rfact_tke * tmask(ji,jj,jk)
280               !                                               ! lower diagonal
281               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) )
282               !                                               ! upper diagonal
283               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) )
284               !                                               ! diagonal
285               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
286               !                                               ! right hand side in en
287               en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
288            END DO
289         END DO
290      END DO
291      !
292      zdiag(:,:,jpk) = 1._wp
293      !
294      ! Set surface condition on zwall_psi (1 at the bottom)
295      zwall_psi(:,:, 1 ) = zwall_psi(:,:,2)
296      zwall_psi(:,:,jpk) = 1._wp
297      !
298      ! Surface boundary condition on tke
299      ! ---------------------------------
300      !
301      SELECT CASE ( nn_bc_surf )
302      !
303      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2)
304      ! First level
305      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3  )
306      zd_lw(:,:,1) = en(:,:,1)
307      zd_up(:,:,1) = 0._wp
308      zdiag(:,:,1) = 1._wp
309      !
310      ! One level below
311      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))   &
312         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   )
313      zd_lw(:,:,2) = 0._wp 
314      zd_up(:,:,2) = 0._wp
315      zdiag(:,:,2) = 1._wp
316      !
317      !
318      CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
319      !
320      ! Dirichlet conditions at k=1
321      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin  )
322      zd_lw(:,:,1) = en(:,:,1)
323      zd_up(:,:,1) = 0._wp
324      zdiag(:,:,1) = 1._wp
325      !
326      ! at k=2, set de/dz=Fw
327      !cbr
328      zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
329      zd_lw(:,:,2) = 0._wp
330      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
331      zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
332          &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
333!!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf)
334      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
335      !
336      !
337      END SELECT
338
339      ! Bottom boundary condition on tke
340      ! --------------------------------
341      !
342      SELECT CASE ( nn_bc_bot )
343      !
344      CASE ( 0 )             ! Dirichlet
345         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
346         !                      ! Balance between the production and the dissipation terms
347         DO jj = 2, jpjm1
348            DO ji = fs_2, fs_jpim1   ! vector opt.
349!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
350!!   With thick deep ocean level thickness, this may be quite large, no ???
351!!   in particular in ocean cavities where top stratification can be large...
352               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
353               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
354               !
355               z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
356               !
357               ! Dirichlet condition applied at:
358               !     Bottom level (ibot)      &      Just above it (ibotm1)   
359               zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp
360               zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
361               zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp
362               en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en
363            END DO
364         END DO
365         !
366         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
367            DO jj = 2, jpjm1
368               DO ji = fs_2, fs_jpim1   ! vector opt.
369                  itop   = mikt(ji,jj)       ! k   top w-point
370                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
371                  !                                                ! mask at the ocean surface points
372                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
373                  !
374 !!gm TO BE VERIFIED !!!
375                  ! Dirichlet condition applied at:
376                  !     top level (itop)         &      Just below it (itopp1)   
377                  zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
378                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
379                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
380                  en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
381               END DO
382            END DO
383         ENDIF
384         !
385      CASE ( 1 )             ! Neumman boundary condition
386         !                     
387         DO jj = 2, jpjm1
388            DO ji = fs_2, fs_jpim1   ! vector opt.
389               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
390               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
391               !
392               z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
393               !
394               ! Bottom level Dirichlet condition:
395               !     Bottom level (ibot)      &      Just above it (ibotm1)   
396               !         Dirichlet            !         Neumann
397               zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag
398               zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
399               zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
400            END DO
401         END DO
402         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
403            DO jj = 2, jpjm1
404               DO ji = fs_2, fs_jpim1   ! vector opt.
405                  itop   = mikt(ji,jj)       ! k   top w-point
406                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
407                  !                                                ! mask at the ocean surface points
408                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
409                  !
410                  ! Bottom level Dirichlet condition:
411                  !     Bottom level (ibot)      &      Just above it (ibotm1)   
412                  !         Dirichlet            !         Neumann
413                  zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
414                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
415                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
416               END DO
417            END DO
418         ENDIF
419         !
420      END SELECT
421
422      ! Matrix inversion (en prescribed at surface and the bottom)
423      ! ----------------------------------------------------------
424      !
425      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
426         DO jj = 2, jpjm1
427            DO ji = fs_2, fs_jpim1    ! vector opt.
428               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
429            END DO
430         END DO
431      END DO
432      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
433         DO jj = 2, jpjm1
434            DO ji = fs_2, fs_jpim1    ! vector opt.
435               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)
436            END DO
437         END DO
438      END DO
439      DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
440         DO jj = 2, jpjm1
441            DO ji = fs_2, fs_jpim1    ! vector opt.
442               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
443            END DO
444         END DO
445      END DO
446      !                                            ! set the minimum value of tke
447      en(:,:,:) = MAX( en(:,:,:), rn_emin )
448
449      !!----------------------------------------!!
450      !!   Solve prognostic equation for psi    !!
451      !!----------------------------------------!!
452
453      ! Set psi to previous time step value
454      !
455      SELECT CASE ( nn_clos )
456      !
457      CASE( 0 )               ! k-kl  (Mellor-Yamada)
458         DO jk = 2, jpkm1
459            DO jj = 2, jpjm1
460               DO ji = fs_2, fs_jpim1   ! vector opt.
461                  psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
462               END DO
463            END DO
464         END DO
465         !
466      CASE( 1 )               ! k-eps
467         DO jk = 2, jpkm1
468            DO jj = 2, jpjm1
469               DO ji = fs_2, fs_jpim1   ! vector opt.
470                  psi(ji,jj,jk)  = eps(ji,jj,jk)
471               END DO
472            END DO
473         END DO
474         !
475      CASE( 2 )               ! k-w
476         DO jk = 2, jpkm1
477            DO jj = 2, jpjm1
478               DO ji = fs_2, fs_jpim1   ! vector opt.
479                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
480               END DO
481            END DO
482         END DO
483         !
484      CASE( 3 )               ! generic
485         DO jk = 2, jpkm1
486            DO jj = 2, jpjm1
487               DO ji = fs_2, fs_jpim1   ! vector opt.
488                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
489               END DO
490            END DO
491         END DO
492         !
493      END SELECT
494      !
495      ! Now gls (output in psi)
496      ! -------------------------------
497      ! Resolution of a tridiagonal linear system by a "methode de chasse"
498      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
499      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
500      ! Warning : after this step, en : right hand side of the matrix
501
502      DO jk = 2, jpkm1
503         DO jj = 2, jpjm1
504            DO ji = fs_2, fs_jpim1   ! vector opt.
505               !
506               ! psi / k
507               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
508               !
509               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable)
510               zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
511               !
512               rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p
513               !
514               ! shear prod. - stratif. destruction
515               prod = rpsi1 * zratio * p_sh2(ji,jj,jk)
516               !
517               ! stratif. destruction
518               buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
519               !
520               ! shear prod. - stratif. destruction
521               diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
522               !
523               zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
524               !
525               zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term
526               zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
527               !                                                       
528               ! building the matrix
529               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
530               !                                               ! lower diagonal
531               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) )
532               !                                               ! upper diagonal
533               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) )
534               !                                               ! diagonal
535               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
536               !                                               ! right hand side in psi
537               psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
538            END DO
539         END DO
540      END DO
541      !
542      zdiag(:,:,jpk) = 1._wp
543
544      ! Surface boundary condition on psi
545      ! ---------------------------------
546      !
547      SELECT CASE ( nn_bc_surf )
548      !
549      CASE ( 0 )             ! Dirichlet boundary conditions
550         !
551         ! Surface value
552         zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic
553         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
554         zd_lw(:,:,1) = psi(:,:,1)
555         zd_up(:,:,1) = 0._wp
556         zdiag(:,:,1) = 1._wp
557         !
558         ! One level below
559         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) )))
560         zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:)
561         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
562         zd_lw(:,:,2) = 0._wp
563         zd_up(:,:,2) = 0._wp
564         zdiag(:,:,2) = 1._wp
565         !
566      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
567         !
568         ! Surface value: Dirichlet
569         zdep    (:,:)   = zhsro(:,:) * rl_sf
570         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
571         zd_lw(:,:,1) = psi(:,:,1)
572         zd_up(:,:,1) = 0._wp
573         zdiag(:,:,1) = 1._wp
574         !
575         ! Neumann condition at k=2
576         zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
577         zd_lw(:,:,2) = 0._wp
578         !
579         ! Set psi vertical flux at the surface:
580         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
581         zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
582         zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
583         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &
584            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.)
585         zflxs(:,:)   = zdep(:,:) * zflxs(:,:)
586         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
587         !
588      END SELECT
589
590      ! Bottom boundary condition on psi
591      ! --------------------------------
592      !
593!!gm should be done for ISF (top boundary cond.)
594!!gm so, totally new staff needed      ===>>> think about that !
595!
596      SELECT CASE ( nn_bc_bot )     ! bottom boundary
597      !
598      CASE ( 0 )             ! Dirichlet
599         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
600         !                      ! Balance between the production and the dissipation terms
601         DO jj = 2, jpjm1
602            DO ji = fs_2, fs_jpim1   ! vector opt.
603               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
604               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
605               zdep(ji,jj) = vkarmn * r_z0_bot
606               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
607               zd_lw(ji,jj,ibot) = 0._wp
608               zd_up(ji,jj,ibot) = 0._wp
609               zdiag(ji,jj,ibot) = 1._wp
610               !
611               ! Just above last level, Dirichlet condition again (GOTM like)
612               zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) )
613               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
614               zd_lw(ji,jj,ibotm1) = 0._wp
615               zd_up(ji,jj,ibotm1) = 0._wp
616               zdiag(ji,jj,ibotm1) = 1._wp
617            END DO
618         END DO
619         !
620      CASE ( 1 )             ! Neumman boundary condition
621         !                     
622         DO jj = 2, jpjm1
623            DO ji = fs_2, fs_jpim1   ! vector opt.
624               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
625               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
626               !
627               ! Bottom level Dirichlet condition:
628               zdep(ji,jj) = vkarmn * r_z0_bot
629               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
630               !
631               zd_lw(ji,jj,ibot) = 0._wp
632               zd_up(ji,jj,ibot) = 0._wp
633               zdiag(ji,jj,ibot) = 1._wp
634               !
635               ! Just above last level: Neumann condition with flux injection
636               zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag
637               zd_up(ji,jj,ibotm1) = 0.
638               !
639               ! Set psi vertical flux at the bottom:
640               zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1)
641               zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
642                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
643               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1)
644            END DO
645         END DO
646         !
647      END SELECT
648
649      ! Matrix inversion
650      ! ----------------
651      !
652      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
653         DO jj = 2, jpjm1
654            DO ji = fs_2, fs_jpim1    ! vector opt.
655               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
656            END DO
657         END DO
658      END DO
659      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
660         DO jj = 2, jpjm1
661            DO ji = fs_2, fs_jpim1    ! vector opt.
662               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)
663            END DO
664         END DO
665      END DO
666      DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
667         DO jj = 2, jpjm1
668            DO ji = fs_2, fs_jpim1    ! vector opt.
669               psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
670            END DO
671         END DO
672      END DO
673
674      ! Set dissipation
675      !----------------
676
677      SELECT CASE ( nn_clos )
678      !
679      CASE( 0 )               ! k-kl  (Mellor-Yamada)
680         DO jk = 1, jpkm1
681            DO jj = 2, jpjm1
682               DO ji = fs_2, fs_jpim1   ! vector opt.
683                  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)
684               END DO
685            END DO
686         END DO
687         !
688      CASE( 1 )               ! k-eps
689         DO jk = 1, jpkm1
690            DO jj = 2, jpjm1
691               DO ji = fs_2, fs_jpim1   ! vector opt.
692                  eps(ji,jj,jk) = psi(ji,jj,jk)
693               END DO
694            END DO
695         END DO
696         !
697      CASE( 2 )               ! k-w
698         DO jk = 1, jpkm1
699            DO jj = 2, jpjm1
700               DO ji = fs_2, fs_jpim1   ! vector opt.
701                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
702               END DO
703            END DO
704         END DO
705         !
706      CASE( 3 )               ! generic
707         zcoef = rc0**( 3._wp  + rpp/rnn )
708         zex1  =      ( 1.5_wp + rmm/rnn )
709         zex2  = -1._wp / rnn
710         DO jk = 1, jpkm1
711            DO jj = 2, jpjm1
712               DO ji = fs_2, fs_jpim1   ! vector opt.
713                  eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
714               END DO
715            END DO
716         END DO
717         !
718      END SELECT
719
720      ! Limit dissipation rate under stable stratification
721      ! --------------------------------------------------
722      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time
723         DO jj = 2, jpjm1
724            DO ji = fs_2, fs_jpim1    ! vector opt.
725               ! limitation
726               eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
727               hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
728               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
729               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
730               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) )
731            END DO
732         END DO
733      END DO 
734
735      !
736      ! Stability function and vertical viscosity and diffusivity
737      ! ---------------------------------------------------------
738      !
739      SELECT CASE ( nn_stab_func )
740      !
741      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
742         DO jk = 2, jpkm1
743            DO jj = 2, jpjm1
744               DO ji = fs_2, fs_jpim1   ! vector opt.
745                  ! zcof =  l²/q²
746                  zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
747                  ! Gh = -N²l²/q²
748                  gh = - rn2(ji,jj,jk) * zcof
749                  gh = MIN( gh, rgh0   )
750                  gh = MAX( gh, rghmin )
751                  ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
752                  sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) )
753                  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)
754                  !
755                  ! Store stability function in zstt and zstm
756                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
757                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
758               END DO
759            END DO
760         END DO
761         !
762      CASE ( 2, 3 )               ! Canuto stability functions
763         DO jk = 2, jpkm1
764            DO jj = 2, jpjm1
765               DO ji = fs_2, fs_jpim1   ! vector opt.
766                  ! zcof =  l²/q²
767                  zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
768                  ! Gh = -N²l²/q²
769                  gh = - rn2(ji,jj,jk) * zcof
770                  gh = MIN( gh, rgh0   )
771                  gh = MAX( gh, rghmin )
772                  gh = gh * rf6
773                  ! Gm =  M²l²/q² Shear number
774                  shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
775                  gm = MAX( shr * zcof , 1.e-10 )
776                  gm = gm * rf6
777                  gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm )
778                  ! Stability functions from Canuto
779                  rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm
780                  sm = (rs0 - rs1*gh + rs2*gm) / rcff
781                  sh = (rs4 - rs5*gh + rs6*gm) / rcff
782                  !
783                  ! Store stability function in zstt and zstm
784                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
785                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
786               END DO
787            END DO
788         END DO
789         !
790      END SELECT
791
792      ! Boundary conditions on stability functions for momentum (Neumann):
793      ! Lines below are useless if GOTM style Dirichlet conditions are used
794
795      zstm(:,:,1) = zstm(:,:,2)
796
797      DO jj = 2, jpjm1
798         DO ji = fs_2, fs_jpim1   ! vector opt.
799            zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
800         END DO
801      END DO
802!!gm should be done for ISF (top boundary cond.)
803!!gm so, totally new staff needed!!gm
804
805      ! Compute diffusivities/viscosities
806      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
807      DO jk = 1, jpk
808         DO jj = 2, jpjm1
809            DO ji = fs_2, fs_jpim1   ! vector opt.
810               zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
811               zavt  = zsqen * zstt(ji,jj,jk)
812               zavm  = zsqen * zstm(ji,jj,jk)
813               p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
814               p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
815            END DO
816         END DO
817      END DO
818      p_avt(:,:,1) = 0._wp
819      !
820      IF(ln_ctl) THEN
821         CALL prt_ctl( tab3d_1=en   , clinfo1=' gls  - e: ', tab3d_2=p_avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
822         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls  - m: ', ovlap=1, kdim=jpk )
823      ENDIF
824      !
825      IF( ln_timing )   CALL timing_stop('zdf_gls')
826      !
827   END SUBROUTINE zdf_gls
828
829
830   SUBROUTINE zdf_gls_init
831      !!----------------------------------------------------------------------
832      !!                  ***  ROUTINE zdf_gls_init  ***
833      !!                     
834      !! ** Purpose :   Initialization of the vertical eddy diffivity and
835      !!              viscosity computed using a GLS turbulent closure scheme
836      !!
837      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
838      !!
839      !! ** input   :   Namlist namzdf_gls
840      !!
841      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
842      !!
843      !!----------------------------------------------------------------------
844      INTEGER ::   jk    ! dummy loop indices
845      INTEGER ::   ios   ! Local integer output status for namelist read
846      REAL(wp)::   zcr   ! local scalar
847      !!
848      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
849         &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
850         &            rn_crban, rn_charn, rn_frac_hs,        &
851         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
852         &            nn_stab_func, nn_clos
853      !!----------------------------------------------------------
854      !
855      IF( ln_timing )   CALL timing_start('zdf_gls_init')
856      !
857      REWIND( numnam_ref )              ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
858      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
859901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in reference namelist', lwp )
860
861      REWIND( numnam_cfg )              ! Namelist namzdf_gls in configuration namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
862      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
863902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp )
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( ln_timing )   CALL timing_stop('zdf_gls_init')
1133      !
1134   END SUBROUTINE zdf_gls_init
1135
1136
1137   SUBROUTINE gls_rst( kt, cdrw )
1138      !!---------------------------------------------------------------------
1139      !!                   ***  ROUTINE ts_rst  ***
1140      !!                     
1141      !! ** Purpose :   Read or write TKE file (en) in restart file
1142      !!
1143      !! ** Method  :   use of IOM library
1144      !!                if the restart does not contain TKE, en is either
1145      !!                set to rn_emin or recomputed (nn_igls/=0)
1146      !!----------------------------------------------------------------------
1147      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
1148      !!
1149      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1150      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1151      !
1152      INTEGER ::   jit, jk   ! dummy loop indices
1153      INTEGER ::   id1, id2, id3, id4
1154      INTEGER ::   ji, jj, ikbu, ikbv
1155      REAL(wp)::   cbx, cby
1156      !!----------------------------------------------------------------------
1157      !
1158      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1159         !                                   ! ---------------
1160         IF( ln_rstart ) THEN                   !* Read the restart file
1161            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. )
1162            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
1163            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
1164            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
1165            !
1166            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
1167               CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
1168               CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k  )
1169               CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k  )
1170               CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n )
1171            ELSE                       
1172               IF(lwp) WRITE(numout,*)
1173               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values'
1174               en    (:,:,:) = rn_emin
1175               hmxl_n(:,:,:) = 0.05_wp
1176               ! avt_k, avm_k already set to the background value in zdf_phy_init
1177            ENDIF
1178         ELSE                                   !* Start from rest
1179            IF(lwp) WRITE(numout,*)
1180            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by 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         !
1186      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1187         !                                   ! -------------------
1188         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1189         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en     ) 
1190         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k  )
1191         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k  )
1192         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
1193         !
1194      ENDIF
1195      !
1196   END SUBROUTINE gls_rst
1197
1198   !!======================================================================
1199END MODULE zdfgls
1200
Note: See TracBrowser for help on using the repository browser.