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

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

source: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/ZDF/zdfgls.F90 @ 13463

Last change on this file since 13463 was 13463, checked in by andmirek, 4 years ago

Ticket #2195:update to trunk 13461

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