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

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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ZDF/zdfgls.F90 @ 11949

Last change on this file since 11949 was 11949, checked in by acc, 4 years ago

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