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/trunk/src/OCE/ZDF – NEMO

source: NEMO/trunk/src/OCE/ZDF/zdfgls.F90 @ 13506

Last change on this file since 13506 was 13506, checked in by smasson, 20 months ago

trunk: avoid unnecessary division by 0 in gls

  • Property svn:keywords set to Id
File size: 57.3 KB
Line 
1MODULE zdfgls
2   !!======================================================================
3   !!                       ***  MODULE  zdfgls  ***
4   !! Ocean physics:  vertical mixing coefficient computed from the gls
5   !!                 turbulent closure parameterization
6   !!======================================================================
7   !! History :  3.0  !  2009-09  (G. Reffray)  Original code
8   !!            3.3  !  2010-10  (C. Bricaud)  Add in the reference
9   !!            4.0  !  2017-04  (G. Madec)  remove CPP keys & avm at t-point only
10   !!             -   !  2017-05  (G. Madec)  add top friction as boundary condition
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   zdf_gls       : update momentum and tracer Kz from a gls scheme
15   !!   zdf_gls_init  : initialization, namelist read, and parameters control
16   !!   gls_rst       : read/write gls restart in ocean restart file
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and active tracers
19   USE dom_oce        ! ocean space and time domain
20   USE domvvl         ! ocean space and time domain : variable volume layer
21   USE zdfdrg  , ONLY : 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      zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
330      zd_lw(:,:,2) = 0._wp
331      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) ))
332      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
333          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
334!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf)
335      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
336      !
337      !
338      END SELECT
339
340      ! Bottom boundary condition on tke
341      ! --------------------------------
342      !
343      SELECT CASE ( nn_bc_bot )
344      !
345      CASE ( 0 )             ! Dirichlet
346         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
347         !                      ! Balance between the production and the dissipation terms
348         DO_2D( 0, 0, 0, 0 )
349!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
350!!   With thick deep ocean level thickness, this may be quite large, no ???
351!!   in particular in ocean cavities where top stratification can be large...
352            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
353            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
354            !
355            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
356            !
357            ! Dirichlet condition applied at:
358            !     Bottom level (ibot)      &      Just above it (ibotm1)   
359            zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp
360            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
361            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp
362            en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en
363         END_2D
364         !
365         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
366            DO_2D( 0, 0, 0, 0 )
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_2D
380         ENDIF
381         !
382      CASE ( 1 )             ! Neumman boundary condition
383         !                     
384         DO_2D( 0, 0, 0, 0 )
385            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
386            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
387            !
388            z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
389            !
390            ! Bottom level Dirichlet condition:
391            !     Bottom level (ibot)      &      Just above it (ibotm1)   
392            !         Dirichlet            !         Neumann
393            zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag
394            zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
395            zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
396         END_2D
397         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
398            DO_2D( 0, 0, 0, 0 )
399               itop   = mikt(ji,jj)       ! k   top w-point
400               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
401               !                                                ! mask at the ocean surface points
402               z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
403               !
404               ! Bottom level Dirichlet condition:
405               !     Bottom level (ibot)      &      Just above it (ibotm1)   
406               !         Dirichlet            !         Neumann
407               zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
408               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
409               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
410            END_2D
411         ENDIF
412         !
413      END SELECT
414
415      ! Matrix inversion (en prescribed at surface and the bottom)
416      ! ----------------------------------------------------------
417      !
418      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
419         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
420      END_3D
421      DO_3D( 0, 0, 0, 0, 2, jpk )                  ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
422         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)
423      END_3D
424      DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
425         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
426      END_3D
427      !                                            ! set the minimum value of tke
428      en(:,:,:) = MAX( en(:,:,:), rn_emin )
429
430      !!----------------------------------------!!
431      !!   Solve prognostic equation for psi    !!
432      !!----------------------------------------!!
433
434      ! Set psi to previous time step value
435      !
436      SELECT CASE ( nn_clos )
437      !
438      CASE( 0 )               ! k-kl  (Mellor-Yamada)
439         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
440            psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
441         END_3D
442         !
443      CASE( 1 )               ! k-eps
444         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
445            psi(ji,jj,jk)  = eps(ji,jj,jk)
446         END_3D
447         !
448      CASE( 2 )               ! k-w
449         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
450            psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
451         END_3D
452         !
453      CASE( 3 )               ! generic
454         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
455            psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
456         END_3D
457         !
458      END SELECT
459      !
460      ! Now gls (output in psi)
461      ! -------------------------------
462      ! Resolution of a tridiagonal linear system by a "methode de chasse"
463      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
464      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
465      ! Warning : after this step, en : right hand side of the matrix
466
467      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
468         !
469         ! psi / k
470         zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
471         !
472         ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable)
473         zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
474         !
475         rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p
476         !
477         ! shear prod. - stratif. destruction
478         prod = rpsi1 * zratio * p_sh2(ji,jj,jk)
479         !
480         ! stratif. destruction
481         buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
482         !
483         ! shear prod. - stratif. destruction
484         diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
485         !
486         zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
487         !
488         zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term
489         zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
490         !                                                       
491         ! building the matrix
492         zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
493         !                                               ! lower diagonal
494         zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) )   &
495            &            / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
496         !                                               ! upper diagonal
497         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) )   &
498            &            / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) )
499         !                                               ! diagonal
500         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk)
501         !                                               ! right hand side in psi
502         psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk)
503      END_3D
504      !
505      zdiag(:,:,jpk) = 1._wp
506
507      ! Surface boundary condition on psi
508      ! ---------------------------------
509      !
510      SELECT CASE ( nn_bc_surf )
511      !
512      CASE ( 0 )             ! Dirichlet boundary conditions
513         !
514         ! Surface value
515         zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic
516         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
517         zd_lw(:,:,1) = psi(:,:,1)
518         zd_up(:,:,1) = 0._wp
519         zdiag(:,:,1) = 1._wp
520         !
521         ! One level below
522         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) )))
523         zdep    (:,:)   = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:)
524         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
525         zd_lw(:,:,2) = 0._wp
526         zd_up(:,:,2) = 0._wp
527         zdiag(:,:,2) = 1._wp
528         !
529      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
530         !
531         ! Surface value: Dirichlet
532         zdep    (:,:)   = zhsro(:,:) * rl_sf
533         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
534         zd_lw(:,:,1) = psi(:,:,1)
535         zd_up(:,:,1) = 0._wp
536         zdiag(:,:,1) = 1._wp
537         !
538         ! Neumann condition at k=2
539         zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
540         zd_lw(:,:,2) = 0._wp
541         !
542         ! Set psi vertical flux at the surface:
543         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope
544         zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf)
545         zflxs(:,:)   = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) &
546            &           *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
547         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &
548            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.)
549         zflxs(:,:)   = zdep(:,:) * zflxs(:,:)
550         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
551         !
552      END SELECT
553
554      ! Bottom boundary condition on psi
555      ! --------------------------------
556      !
557!!gm should be done for ISF (top boundary cond.)
558!!gm so, totally new staff needed      ===>>> think about that !
559!
560      SELECT CASE ( nn_bc_bot )     ! bottom boundary
561      !
562      CASE ( 0 )             ! Dirichlet
563         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
564         !                      ! Balance between the production and the dissipation terms
565         DO_2D( 0, 0, 0, 0 )
566            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
567            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
568            zdep(ji,jj) = vkarmn * r_z0_bot
569            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
570            zd_lw(ji,jj,ibot) = 0._wp
571            zd_up(ji,jj,ibot) = 0._wp
572            zdiag(ji,jj,ibot) = 1._wp
573            !
574            ! Just above last level, Dirichlet condition again (GOTM like)
575            zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) )
576            psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
577            zd_lw(ji,jj,ibotm1) = 0._wp
578            zd_up(ji,jj,ibotm1) = 0._wp
579            zdiag(ji,jj,ibotm1) = 1._wp
580         END_2D
581         !
582      CASE ( 1 )             ! Neumman boundary condition
583         !                     
584         DO_2D( 0, 0, 0, 0 )
585            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
586            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
587            !
588            ! Bottom level Dirichlet condition:
589            zdep(ji,jj) = vkarmn * r_z0_bot
590            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
591            !
592            zd_lw(ji,jj,ibot) = 0._wp
593            zd_up(ji,jj,ibot) = 0._wp
594            zdiag(ji,jj,ibot) = 1._wp
595            !
596            ! Just above last level: Neumann condition with flux injection
597            zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag
598            zd_up(ji,jj,ibotm1) = 0.
599            !
600            ! Set psi vertical flux at the bottom:
601            zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm)
602            zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
603               &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
604            psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm)
605         END_2D
606         !
607      END SELECT
608
609      ! Matrix inversion
610      ! ----------------
611      !
612      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
613         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
614      END_3D
615      DO_3D( 0, 0, 0, 0, 2, jpk )                  ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
616         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)
617      END_3D
618      DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
619         psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
620      END_3D
621
622      ! Set dissipation
623      !----------------
624
625      SELECT CASE ( nn_clos )
626      !
627      CASE( 0 )               ! k-kl  (Mellor-Yamada)
628         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
629            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)
630         END_3D
631         !
632      CASE( 1 )               ! k-eps
633         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
634            eps(ji,jj,jk) = psi(ji,jj,jk)
635         END_3D
636         !
637      CASE( 2 )               ! k-w
638         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
639            eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
640         END_3D
641         !
642      CASE( 3 )               ! generic
643         zcoef = rc0**( 3._wp  + rpp/rnn )
644         zex1  =      ( 1.5_wp + rmm/rnn )
645         zex2  = -1._wp / rnn
646         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
647            eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
648         END_3D
649         !
650      END SELECT
651
652      ! Limit dissipation rate under stable stratification
653      ! --------------------------------------------------
654      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
655         ! limitation
656         eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
657         hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
658         ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
659         zrn2 = MAX( rn2(ji,jj,jk), rsmall )
660         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) )
661      END_3D
662
663      !
664      ! Stability function and vertical viscosity and diffusivity
665      ! ---------------------------------------------------------
666      !
667      SELECT CASE ( nn_stab_func )
668      !
669      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
670         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
671            ! zcof =  l²/q²
672            zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
673            ! Gh = -N²l²/q²
674            gh = - rn2(ji,jj,jk) * zcof
675            gh = MIN( gh, rgh0   )
676            gh = MAX( gh, rghmin )
677            ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
678            sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) )
679            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)
680            !
681            ! Store stability function in zstt and zstm
682            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
683            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
684         END_3D
685         !
686      CASE ( 2, 3 )               ! Canuto stability functions
687         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
688            ! zcof =  l²/q²
689            zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
690            ! Gh = -N²l²/q²
691            gh = - rn2(ji,jj,jk) * zcof
692            gh = MIN( gh, rgh0   )
693            gh = MAX( gh, rghmin )
694            gh = gh * rf6
695            ! Gm =  M²l²/q² Shear number
696            shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
697            gm = MAX( shr * zcof , 1.e-10 )
698            gm = gm * rf6
699            gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm )
700            ! Stability functions from Canuto
701            rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm
702            sm = (rs0 - rs1*gh + rs2*gm) / rcff
703            sh = (rs4 - rs5*gh + rs6*gm) / rcff
704            !
705            ! Store stability function in zstt and zstm
706            zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
707            zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
708         END_3D
709         !
710      END SELECT
711
712      ! Boundary conditions on stability functions for momentum (Neumann):
713      ! Lines below are useless if GOTM style Dirichlet conditions are used
714
715      zstm(:,:,1) = zstm(:,:,2)
716
717      ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
718      zstm(:,:,jpk) = 0. 
719      DO_2D( 0, 0, 0, 0 )             ! update bottom with good values
720         zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
721      END_2D
722
723      zstt(:,:,  1) = wmask(:,:,  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
724      zstt(:,:,jpk) = wmask(:,:,jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
725
726!!gm should be done for ISF (top boundary cond.)
727!!gm so, totally new staff needed!!gm
728
729      ! Compute diffusivities/viscosities
730      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
731      !  -> yes BUT p_avm(:,:1) and p_avm(:,:jpk) are used when we compute zd_lw(:,:2) and zd_up(:,:jpkm1). These values are
732      !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
733      !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
734      DO_3D( 0, 0, 0, 0, 1, jpk )
735         zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
736         zavt  = zsqen * zstt(ji,jj,jk)
737         zavm  = zsqen * zstm(ji,jj,jk)
738         p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
739         p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
740      END_3D
741      p_avt(:,:,1) = 0._wp
742      !
743      IF(sn_cfctl%l_prtctl) THEN
744         CALL prt_ctl( tab3d_1=en   , clinfo1=' gls  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
745         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls  - m: ', kdim=jpk )
746      ENDIF
747      !
748   END SUBROUTINE zdf_gls
749
750
751   SUBROUTINE zdf_gls_init
752      !!----------------------------------------------------------------------
753      !!                  ***  ROUTINE zdf_gls_init  ***
754      !!                     
755      !! ** Purpose :   Initialization of the vertical eddy diffivity and
756      !!              viscosity computed using a GLS turbulent closure scheme
757      !!
758      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
759      !!
760      !! ** input   :   Namlist namzdf_gls
761      !!
762      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
763      !!
764      !!----------------------------------------------------------------------
765      INTEGER ::   jk    ! dummy loop indices
766      INTEGER ::   ios   ! Local integer output status for namelist read
767      REAL(wp)::   zcr   ! local scalar
768      !!
769      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       &
770         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   &
771         &            rn_crban, rn_charn, rn_frac_hs,              &
772         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, &
773         &            nn_stab_func, nn_clos
774      !!----------------------------------------------------------
775      !
776      READ  ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
777901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_gls in reference namelist' )
778
779      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
780902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist' )
781      IF(lwm) WRITE ( numond, namzdf_gls )
782
783      IF(lwp) THEN                     !* Control print
784         WRITE(numout,*)
785         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
786         WRITE(numout,*) '~~~~~~~~~~~~'
787         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
788         WRITE(numout,*) '      minimum value of en                           rn_emin        = ', rn_emin
789         WRITE(numout,*) '      minimum value of eps                          rn_epsmin      = ', rn_epsmin
790         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim  = ', ln_length_lim
791         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp   = ', rn_clim_galp
792         WRITE(numout,*) '      TKE Surface boundary condition                nn_bc_surf     = ', nn_bc_surf
793         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
794         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
795         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
796         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
797         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
798         WRITE(numout,*) '      surface wave breaking under ice               nn_z0_ice      = ', nn_z0_ice
799         SELECT CASE( nn_z0_ice )
800         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on surface wave breaking'
801         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )'
802         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-fr_i(:,:)'
803         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )'
804         CASE DEFAULT
805            CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3')
806         END SELECT
807         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs
808         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
809         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
810         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
811         WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri
812         WRITE(numout,*)
813         WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
814         WRITE(numout,*) '      top    ocean cavity roughness (m)             rn_z0(_top)   = ', r_z0_top
815         WRITE(numout,*) '      Bottom seafloor     roughness (m)             rn_z0(_bot)   = ', r_z0_bot
816         WRITE(numout,*)
817      ENDIF
818
819      !                                !* allocate GLS arrays
820      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
821
822      !                                !* Check of some namelist values
823      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' )
824      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' )
825      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' )
826      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' )
827      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' )
828      IF( nn_clos      < 0 .OR. nn_clos      > 3 )              CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
829
830      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
831      !
832      CASE( 0 )                              ! k-kl  (Mellor-Yamada)
833         !
834         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
835         IF(lwp) WRITE(numout,*)
836         rpp     = 0._wp
837         rmm     = 1._wp
838         rnn     = 1._wp
839         rsc_tke = 1.96_wp
840         rsc_psi = 1.96_wp
841         rpsi1   = 0.9_wp
842         rpsi3p  = 1._wp
843         rpsi2   = 0.5_wp
844         !
845         SELECT CASE ( nn_stab_func )
846         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
847         CASE( 2 )      ;   rpsi3m = 2.62_wp       ! Canuto A stability functions
848         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
849         END SELECT
850         !
851      CASE( 1 )                              ! k-eps
852         !
853         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen'
854         IF(lwp) WRITE(numout,*)
855         rpp     =  3._wp
856         rmm     =  1.5_wp
857         rnn     = -1._wp
858         rsc_tke =  1._wp
859         rsc_psi =  1.2_wp  ! Schmidt number for psi
860         rpsi1   =  1.44_wp
861         rpsi3p  =  1._wp
862         rpsi2   =  1.92_wp
863         !
864         SELECT CASE ( nn_stab_func )
865         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
866         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
867         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
868         END SELECT
869         !
870      CASE( 2 )                              ! k-omega
871         !
872         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen'
873         IF(lwp) WRITE(numout,*)
874         rpp     = -1._wp
875         rmm     =  0.5_wp
876         rnn     = -1._wp
877         rsc_tke =  2._wp
878         rsc_psi =  2._wp
879         rpsi1   =  0.555_wp
880         rpsi3p  =  1._wp
881         rpsi2   =  0.833_wp
882         !
883         SELECT CASE ( nn_stab_func )
884         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
885         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
886         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
887         END SELECT
888         !
889      CASE( 3 )                              ! generic
890         !
891         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen'
892         IF(lwp) WRITE(numout,*)
893         rpp     = 2._wp
894         rmm     = 1._wp
895         rnn     = -0.67_wp
896         rsc_tke = 0.8_wp
897         rsc_psi = 1.07_wp
898         rpsi1   = 1._wp
899         rpsi3p  = 1._wp
900         rpsi2   = 1.22_wp
901         !
902         SELECT CASE ( nn_stab_func )
903         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
904         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
905         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
906         END SELECT
907         !
908      END SELECT
909
910      !
911      SELECT CASE ( nn_stab_func )     !* set the parameters of the stability functions
912      !
913      CASE ( 0 )                             ! Galperin stability functions
914         !
915         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin'
916         rc2     =  0._wp
917         rc3     =  0._wp
918         rc_diff =  1._wp
919         rc0     =  0.5544_wp
920         rcm_sf  =  0.9884_wp
921         rghmin  = -0.28_wp
922         rgh0    =  0.0233_wp
923         rghcri  =  0.02_wp
924         !
925      CASE ( 1 )                             ! Kantha-Clayson stability functions
926         !
927         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson'
928         rc2     =  0.7_wp
929         rc3     =  0.2_wp
930         rc_diff =  1._wp
931         rc0     =  0.5544_wp
932         rcm_sf  =  0.9884_wp
933         rghmin  = -0.28_wp
934         rgh0    =  0.0233_wp
935         rghcri  =  0.02_wp
936         !
937      CASE ( 2 )                             ! Canuto A stability functions
938         !
939         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A'
940         rs0 = 1.5_wp * rl1 * rl5*rl5
941         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
942         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
943         rs4 = 2._wp * rl5
944         rs5 = 2._wp * rl4
945         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
946            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
947         rd0 = 3._wp * rl5*rl5
948         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
949         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
950         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
951         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
952         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
953         rc0 = 0.5268_wp
954         rf6 = 8._wp / (rc0**6._wp)
955         rc_diff = SQRT(2._wp) / (rc0**3._wp)
956         rcm_sf  =  0.7310_wp
957         rghmin  = -0.28_wp
958         rgh0    =  0.0329_wp
959         rghcri  =  0.03_wp
960         !
961      CASE ( 3 )                             ! Canuto B stability functions
962         !
963         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B'
964         rs0 = 1.5_wp * rm1 * rm5*rm5
965         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
966         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
967         rs4 = 2._wp * rm5
968         rs5 = 2._wp * rm4
969         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)
970         rd0 = 3._wp * rm5*rm5
971         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
972         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
973         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
974         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
975         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
976         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
977         rf6 = 8._wp / ( rc0**6._wp )
978         rc_diff = SQRT(2._wp)/(rc0**3.)
979         rcm_sf  =  0.7470_wp
980         rghmin  = -0.28_wp
981         rgh0    =  0.0444_wp
982         rghcri  =  0.0414_wp
983         !
984      END SELECT
985   
986      !                                !* Set Schmidt number for psi diffusion in the wave breaking case
987      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
988      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
989      IF( ln_sigpsi ) THEN
990         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
991         ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
992         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
993         ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn
994         rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.))
995      ELSE
996         rsc_psi0 = rsc_psi
997      ENDIF
998 
999      !                                !* Shear free turbulence parameters
1000      !
1001      ra_sf  = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) &
1002               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
1003
1004      IF ( rn_crban==0._wp ) THEN
1005         rl_sf = vkarmn
1006      ELSE
1007         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp) * rsc_tke        &
1008            &                                            + 12._wp*rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm)   &
1009            &                                                     *SQRT(rsc_tke*(rsc_tke                 &
1010            &                                                        + 24._wp*rsc_psi0*rpsi2)) )         &
1011            &                                              /(12._wp*rnn**2.)                             )
1012      ENDIF
1013
1014      !
1015      IF(lwp) THEN                     !* Control print
1016         WRITE(numout,*)
1017         WRITE(numout,*) '   Limit values :'
1018         WRITE(numout,*) '      Parameter  m = ', rmm
1019         WRITE(numout,*) '      Parameter  n = ', rnn
1020         WRITE(numout,*) '      Parameter  p = ', rpp
1021         WRITE(numout,*) '      rpsi1    = ', rpsi1
1022         WRITE(numout,*) '      rpsi2    = ', rpsi2
1023         WRITE(numout,*) '      rpsi3m   = ', rpsi3m
1024         WRITE(numout,*) '      rpsi3p   = ', rpsi3p
1025         WRITE(numout,*) '      rsc_tke  = ', rsc_tke
1026         WRITE(numout,*) '      rsc_psi  = ', rsc_psi
1027         WRITE(numout,*) '      rsc_psi0 = ', rsc_psi0
1028         WRITE(numout,*) '      rc0      = ', rc0
1029         WRITE(numout,*)
1030         WRITE(numout,*) '   Shear free turbulence parameters:'
1031         WRITE(numout,*) '      rcm_sf   = ', rcm_sf
1032         WRITE(numout,*) '      ra_sf    = ', ra_sf
1033         WRITE(numout,*) '      rl_sf    = ', rl_sf
1034      ENDIF
1035
1036      !                                !* Constants initialization
1037      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
1038      rc03  = rc02 * rc0
1039      rc04  = rc03 * rc0
1040      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
1041      rsbc_tke2 = rn_Dt * rn_crban / rl_sf                                 ! Neumann + Wave breaking
1042      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
1043      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer
1044      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
1045      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
1046      rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi
1047      rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
1048      !
1049      rfact_tke = -0.5_wp / rsc_tke * rn_Dt                                ! Cst used for the Diffusion term of tke
1050      rfact_psi = -0.5_wp / rsc_psi * rn_Dt                                ! Cst used for the Diffusion term of tke
1051      !
1052      !                                !* Wall proximity function
1053!!gm tmask or wmask ????
1054      zwall(:,:,:) = 1._wp * tmask(:,:,:)
1055
1056      !                                !* read or initialize all required files 
1057      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxl_n)
1058      !
1059      IF( lwxios ) THEN
1060         CALL iom_set_rstw_var_active('en')
1061         CALL iom_set_rstw_var_active('avt_k')
1062         CALL iom_set_rstw_var_active('avm_k')
1063         CALL iom_set_rstw_var_active('hmxl_n')
1064      ENDIF
1065      !
1066   END SUBROUTINE zdf_gls_init
1067
1068
1069   SUBROUTINE gls_rst( kt, cdrw )
1070      !!---------------------------------------------------------------------
1071      !!                   ***  ROUTINE gls_rst  ***
1072      !!                     
1073      !! ** Purpose :   Read or write TKE file (en) in restart file
1074      !!
1075      !! ** Method  :   use of IOM library
1076      !!                if the restart does not contain TKE, en is either
1077      !!                set to rn_emin or recomputed (nn_igls/=0)
1078      !!----------------------------------------------------------------------
1079      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
1080      !!
1081      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1082      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1083      !
1084      INTEGER ::   jit, jk   ! dummy loop indices
1085      INTEGER ::   id1, id2, id3, id4
1086      INTEGER ::   ji, jj, ikbu, ikbv
1087      REAL(wp)::   cbx, cby
1088      !!----------------------------------------------------------------------
1089      !
1090      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1091         !                                   ! ---------------
1092         IF( ln_rstart ) THEN                   !* Read the restart file
1093            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. )
1094            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
1095            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
1096            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
1097            !
1098            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
1099               CALL iom_get( numror, jpdom_auto, 'en'    , en    , ldxios = lrxios )
1100               CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k , ldxios = lrxios )
1101               CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k , ldxios = lrxios )
1102               CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n, ldxios = lrxios )
1103            ELSE                       
1104               IF(lwp) WRITE(numout,*)
1105               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values'
1106               en    (:,:,:) = rn_emin
1107               hmxl_n(:,:,:) = 0.05_wp
1108               ! avt_k, avm_k already set to the background value in zdf_phy_init
1109            ENDIF
1110         ELSE                                   !* Start from rest
1111            IF(lwp) WRITE(numout,*)
1112            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values'
1113            en    (:,:,:) = rn_emin
1114            hmxl_n(:,:,:) = 0.05_wp
1115            ! avt_k, avm_k already set to the background value in zdf_phy_init
1116         ENDIF
1117         !
1118      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1119         !                                   ! -------------------
1120         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1121         IF( lwxios ) CALL iom_swap(      cwxios_context         )
1122         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en    , ldxios = lwxios )
1123         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k , ldxios = lwxios )
1124         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k , ldxios = lwxios )
1125         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n, ldxios = lwxios )
1126         IF( lwxios ) CALL iom_swap(      cxios_context          )
1127         !
1128      ENDIF
1129      !
1130   END SUBROUTINE gls_rst
1131
1132   !!======================================================================
1133END MODULE zdfgls
Note: See TracBrowser for help on using the repository browser.