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/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90 @ 13553

Last change on this file since 13553 was 13553, checked in by hadcv, 3 years ago

Merge in trunk up to [13550]

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