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 @ 12236

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

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