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/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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