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 branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago

dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation

  • Property svn:keywords set to Id
File size: 59.6 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   !!----------------------------------------------------------------------
10#if defined key_zdfgls   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_zdfgls'                 Generic Length Scale vertical physics
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 zdf_oce        ! ocean vertical physics
22   USE sbc_oce        ! surface boundary condition: ocean
23   USE phycst         ! physical constants
24   USE zdfmxl         ! mixed layer
25   USE restart        ! only for lrst_oce
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27   USE lib_mpp        ! MPP manager
28   USE prtctl         ! Print control
29   USE in_out_manager ! I/O manager
30   USE iom            ! I/O manager library
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   zdf_gls        ! routine called in step module
36   PUBLIC   zdf_gls_init   ! routine called in opa module
37   PUBLIC   gls_rst        ! routine called in step module
38
39   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfgls = .TRUE.   !: TKE vertical mixing flag
40   !
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en      !: now turbulent kinetic energy
42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points
46
47   !                                         !!! ** Namelist  namzdf_gls  **
48   LOGICAL  ::   ln_crban      = .FALSE.      ! =T use Craig and Banner scheme
49   LOGICAL  ::   ln_length_lim = .FALSE.      ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
50   LOGICAL  ::   ln_sigpsi     = .FALSE.      ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
51   INTEGER  ::   nn_tkebc_surf = 0            ! TKE surface boundary condition (=0/1)
52   INTEGER  ::   nn_tkebc_bot  = 0            ! TKE bottom boundary condition (=0/1)
53   INTEGER  ::   nn_psibc_surf = 0            ! PSI surface boundary condition (=0/1)
54   INTEGER  ::   nn_psibc_bot  = 0            ! PSI bottom boundary condition (=0/1)
55   INTEGER  ::   nn_stab_func  = 0            ! stability functions G88, KC or Canuto (=0/1/2)
56   INTEGER  ::   nn_clos       = 0            ! closure 0/1/2/3 MY82/k-eps/k-w/gen
57   REAL(wp) ::   rn_clim_galp  = 0.53_wp      ! Holt 2008 value for k-eps: 0.267
58   REAL(wp) ::   rn_epsmin     = 1.e-12_wp    ! minimum value of dissipation (m2/s3)
59   REAL(wp) ::   rn_emin       = 1.e-6_wp     ! minimum value of TKE (m2/s2)
60   REAL(wp) ::   rn_charn      = 2.e+5_wp     ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
61   REAL(wp) ::   rn_crban      = 100._wp      ! Craig and Banner constant for surface breaking waves mixing
62
63   REAL(wp) ::   hsro          =  0.003_wp    ! Minimum surface roughness
64   REAL(wp) ::   hbro          =  0.003_wp    ! Bottom roughness (m)
65   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
66   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
67   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn   
68   REAL(wp) ::   rghmin        = -0.28_wp
69   REAL(wp) ::   rgh0          =  0.0329_wp
70   REAL(wp) ::   rghcri        =  0.03_wp
71   REAL(wp) ::   ra1           =  0.92_wp
72   REAL(wp) ::   ra2           =  0.74_wp
73   REAL(wp) ::   rb1           = 16.60_wp
74   REAL(wp) ::   rb2           = 10.10_wp         
75   REAL(wp) ::   re2           =  1.33_wp         
76   REAL(wp) ::   rl1           =  0.107_wp
77   REAL(wp) ::   rl2           =  0.0032_wp
78   REAL(wp) ::   rl3           =  0.0864_wp
79   REAL(wp) ::   rl4           =  0.12_wp
80   REAL(wp) ::   rl5           = 11.9_wp
81   REAL(wp) ::   rl6           =  0.4_wp
82   REAL(wp) ::   rl7           =  0.0_wp
83   REAL(wp) ::   rl8           =  0.48_wp
84   REAL(wp) ::   rm1           =  0.127_wp
85   REAL(wp) ::   rm2           =  0.00336_wp
86   REAL(wp) ::   rm3           =  0.0906_wp
87   REAL(wp) ::   rm4           =  0.101_wp
88   REAL(wp) ::   rm5           = 11.2_wp
89   REAL(wp) ::   rm6           =  0.4_wp
90   REAL(wp) ::   rm7           =  0.0_wp
91   REAL(wp) ::   rm8           =  0.318_wp
92   
93   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
94   REAL(wp) ::   rc03_sqrt2_galp                                  !     -           -           -        -
95   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke       !     -           -           -        -
96   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi       !     -           -           -        -
97   REAL(wp) ::   rsbc_mb  , rsbc_std , rsbc_zs                    !     -           -           -        -
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   !! * Substitutions
105#  include "domzgr_substitute.h90"
106#  include "vectopt_loop_substitute.h90"
107   !!----------------------------------------------------------------------
108   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
109   !! $Id$
110   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
111   !!----------------------------------------------------------------------
112CONTAINS
113
114   INTEGER FUNCTION zdf_gls_alloc()
115      !!----------------------------------------------------------------------
116      !!                ***  FUNCTION zdf_gls_alloc  ***
117      !!----------------------------------------------------------------------
118      ALLOCATE( en(jpi,jpj,jpk),  mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
119         &      ustars2(jpi,jpj), ustarb2(jpi,jpj)                      , STAT= zdf_gls_alloc )
120         !
121      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
122      IF( zdf_gls_alloc /= 0 )   CALL ctl_warn('zdf_gls_alloc: failed to allocate arrays')
123   END FUNCTION zdf_gls_alloc
124
125
126   SUBROUTINE zdf_gls( kt )
127      !!----------------------------------------------------------------------
128      !!                   ***  ROUTINE zdf_gls  ***
129      !!
130      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
131      !!              coefficients using the GLS turbulent closure scheme.
132      !!----------------------------------------------------------------------
133      USE oce,     z_elem_a  =>   ua   ! use ua as workspace
134      USE oce,     z_elem_b  =>   va   ! use va as workspace
135      USE oce,     z_elem_c  =>   ta   ! use ta as workspace
136      USE oce,     psi       =>   sa   ! use sa as workspace
137      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
138      USE wrk_nemo, ONLY: zdep  => wrk_2d_1
139      USE wrk_nemo, ONLY: zflxs => wrk_2d_2     ! Turbulence fluxed induced by internal waves
140      USE wrk_nemo, ONLY: zhsro => wrk_2d_3     ! Surface roughness (surface waves)
141      USE wrk_nemo, ONLY: eb        => wrk_3d_1   ! tke at time before
142      USE wrk_nemo, ONLY: mxlb      => wrk_3d_2   ! mixing length at time before
143      USE wrk_nemo, ONLY: shear     => wrk_3d_3   ! vertical shear
144      USE wrk_nemo, ONLY: eps       => wrk_3d_4   ! dissipation rate
145      USE wrk_nemo, ONLY: zwall_psi => wrk_3d_5   ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T)
146      !
147      INTEGER, INTENT(in) ::   kt ! ocean time step
148      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
149      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars
150      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -
151      REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      -
152      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
153      REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      -
154      !!--------------------------------------------------------------------
155
156      IF(  wrk_in_use(2, 1,2,3)  .OR.  wrk_in_use(3, 1,2,3,4,5)  ) THEN
157         CALL ctl_stop('zdf_gls: requested workspace arrays unavailable.')   ;   RETURN
158      END IF
159
160      ! Preliminary computing
161
162      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
163
164      ! Compute surface and bottom friction at T-points
165!CDIR NOVERRCHK
166      DO jj = 2, jpjm1
167!CDIR NOVERRCHK
168         DO ji = fs_2, fs_jpim1   ! vector opt.
169            !
170            ! surface friction
171            ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1)
172            !
173            ! bottom friction (explicit before friction)
174            ! Note that we chose here not to bound the friction as in dynbfr)
175            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &
176               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )
177            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &
178               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )
179            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
180         END DO
181      END DO 
182
183      ! In case of breaking surface waves mixing,
184      ! Compute surface roughness length according to Charnock formula:
185      IF( ln_crban ) THEN   ;   zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro)
186      ELSE                  ;   zhsro(:,:) = hsro
187      ENDIF
188
189      ! Compute shear and dissipation rate
190      DO jk = 2, jpkm1
191         DO jj = 2, jpjm1
192            DO ji = fs_2, fs_jpim1   ! vector opt.
193               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
194                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &
195                  &                            / (  fse3uw_n(ji,jj,jk)               &
196                  &                            *    fse3uw_b(ji,jj,jk) )
197               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
198                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
199                  &                            / (  fse3vw_n(ji,jj,jk)               &
200                  &                            *    fse3vw_b(ji,jj,jk) )
201               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
202            END DO
203         END DO
204      END DO
205      !
206      ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
207      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
208
209      ! Save tke at before time step
210      eb  (:,:,:) = en  (:,:,:)
211      mxlb(:,:,:) = mxln(:,:,:)
212
213      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
214         DO jk = 2, jpkm1
215            DO jj = 2, jpjm1 
216               DO ji = fs_2, fs_jpim1   ! vector opt.
217                  zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)
218                  zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbkt(ji,jj)+1) )
219                  zcoef = ( zup / MAX( zdown, rsmall ) )
220                  zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk)
221               END DO
222            END DO
223         END DO
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      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
237      ! Warning : after this step, en : right hand side of the matrix
238
239      DO jk = 2, jpkm1
240         DO jj = 2, jpjm1
241            DO ji = fs_2, fs_jpim1   ! vector opt.
242               !
243               ! shear prod. at w-point weightened by mask
244               shear(ji,jj,jk) =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
245                  &             + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
246               !
247               ! stratif. destruction
248               buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)
249               !
250               ! shear prod. - stratif. destruction
251               diss = eps(ji,jj,jk)
252               !
253               dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
254               !
255               zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk)          ! production term
256               zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
257               !
258               ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
259               ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style)
260               ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
261               ! Otherwise, this should be rsc_psi/rsc_psi0
262               IF( ln_sigpsi ) THEN
263                  zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
264                  zwall_psi(ji,jj,jk) = rsc_psi / (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
265               ELSE
266                  zwall_psi(ji,jj,jk) = 1._wp
267               ENDIF
268               !
269               ! building the matrix
270               zcof = rfact_tke * tmask(ji,jj,jk)
271               !
272               ! lower diagonal
273               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
274                  &                      / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
275               !
276               ! upper diagonal
277               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
278                  &                      / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
279               !
280               ! diagonal
281               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
282                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
283               !
284               ! right hand side in en
285               en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
286            END DO
287         END DO
288      END DO
289      !
290      z_elem_b(:,:,jpk) = 1._wp
291      !
292      ! Set surface condition on zwall_psi (1 at the bottom)
293      IF( ln_sigpsi ) THEN
294         zcoef = rsc_psi / rsc_psi0
295         DO jj = 2, jpjm1
296            DO ji = fs_2, fs_jpim1   ! vector opt.
297               zwall_psi(ji,jj,1) = zcoef
298            END DO
299         END DO
300      ENDIF
301
302      ! Surface boundary condition on tke
303      ! ---------------------------------
304      !
305      SELECT CASE ( nn_tkebc_surf )
306      !
307      CASE ( 0 )             ! Dirichlet case
308         !
309         IF (ln_crban) THEN     ! Wave induced mixing case
310            !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
311            !                      ! balance between the production and the dissipation terms including the wave effect
312            en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
313            z_elem_a(:,:,1) = en(:,:,1)
314            z_elem_c(:,:,1) = 0._wp
315            z_elem_b(:,:,1) = 1._wp
316            !
317            ! one level below
318            en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
319            z_elem_a(:,:,2) = 0._wp
320            z_elem_c(:,:,2) = 0._wp
321            z_elem_b(:,:,2) = 1._wp
322            !
323         ELSE                   ! No wave induced mixing case
324            !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs
325            !                      ! balance between the production and the dissipation terms
326            en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
327            z_elem_a(:,:,1) = en(:,:,1) 
328            z_elem_c(:,:,1) = 0._wp
329            z_elem_b(:,:,1) = 1._wp
330            !
331            ! one level below
332            en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin )
333            z_elem_a(:,:,2) = 0._wp
334            z_elem_c(:,:,2) = 0._wp
335            z_elem_b(:,:,2) = 1._wp
336            !
337         ENDIF
338         !
339      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
340         !
341         IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw
342            !
343            ! Dirichlet conditions at k=1 (Cosmetic)
344            en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
345            z_elem_a(:,:,1) = en(:,:,1)
346            z_elem_c(:,:,1) = 0._wp
347            z_elem_b(:,:,1) = 1._wp
348            ! at k=2, set de/dz=Fw
349            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
350            z_elem_a(:,:,2) = 0._wp       
351            zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)
352            en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
353            !
354         ELSE                   ! No wave induced mixing case: d(e)/dz=0.
355            !
356            ! Dirichlet conditions at k=1 (Cosmetic)
357            en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
358            z_elem_a(:,:,1) = en(:,:,1)
359            z_elem_c(:,:,1) = 0._wp
360            z_elem_b(:,:,1) = 1._wp
361            ! at k=2 set de/dz=0.:
362            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2)  ! Remove z_elem_a from z_elem_b
363            z_elem_a(:,:,2) = 0._wp
364            !
365         ENDIF
366         !
367      END SELECT
368
369      ! Bottom boundary condition on tke
370      ! --------------------------------
371      !
372      SELECT CASE ( nn_tkebc_bot )
373      !
374      CASE ( 0 )             ! Dirichlet
375         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
376         !                      ! Balance between the production and the dissipation terms
377!CDIR NOVERRCHK
378         DO jj = 2, jpjm1
379!CDIR NOVERRCHK
380            DO ji = fs_2, fs_jpim1   ! vector opt.
381               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
382               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
383               !
384               ! Bottom level Dirichlet condition:
385               z_elem_a(ji,jj,ibot  ) = 0._wp
386               z_elem_c(ji,jj,ibot  ) = 0._wp
387               z_elem_b(ji,jj,ibot  ) = 1._wp
388               en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
389               !
390               ! Just above last level, Dirichlet condition again
391               z_elem_a(ji,jj,ibotm1) = 0._wp
392               z_elem_c(ji,jj,ibotm1) = 0._wp
393               z_elem_b(ji,jj,ibotm1) = 1._wp
394               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
395            END DO
396         END DO
397         !
398      CASE ( 1 )             ! Neumman boundary condition
399         !                     
400!CDIR NOVERRCHK
401         DO jj = 2, jpjm1
402!CDIR NOVERRCHK
403            DO ji = fs_2, fs_jpim1   ! vector opt.
404               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
405               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
406               !
407               ! Bottom level Dirichlet condition:
408               z_elem_a(ji,jj,ibot) = 0._wp
409               z_elem_c(ji,jj,ibot) = 0._wp
410               z_elem_b(ji,jj,ibot) = 1._wp
411               en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
412               !
413               ! Just above last level: Neumann condition
414               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)   ! Remove z_elem_c from z_elem_b
415               z_elem_c(ji,jj,ibotm1) = 0._wp
416            END DO
417         END DO
418         !
419      END SELECT
420
421      ! Matrix inversion (en prescribed at surface and the bottom)
422      ! ----------------------------------------------------------
423      !
424      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
425         DO jj = 2, jpjm1
426            DO ji = fs_2, fs_jpim1    ! vector opt.
427               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
428            END DO
429         END DO
430      END DO
431      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
432         DO jj = 2, jpjm1
433            DO ji = fs_2, fs_jpim1    ! vector opt.
434               z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
435            END DO
436         END DO
437      END DO
438      DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
439         DO jj = 2, jpjm1
440            DO ji = fs_2, fs_jpim1    ! vector opt.
441               en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
442            END DO
443         END DO
444      END DO
445      !                                            ! set the minimum value of tke
446      en(:,:,:) = MAX( en(:,:,:), rn_emin )
447     
448      !!----------------------------------------!!
449      !!   Solve prognostic equation for psi    !!
450      !!----------------------------------------!!
451
452      ! Set psi to previous time step value
453      !
454      SELECT CASE ( nn_clos )
455      !
456      CASE( 0 )               ! k-kl  (Mellor-Yamada)
457         DO jk = 2, jpkm1
458            DO jj = 2, jpjm1
459               DO ji = fs_2, fs_jpim1   ! vector opt.
460                  psi(ji,jj,jk)  = en(ji,jj,jk) * mxln(ji,jj,jk)
461               END DO
462            END DO
463         END DO
464         !
465      CASE( 1 )               ! k-eps
466         DO jk = 2, jpkm1
467            DO jj = 2, jpjm1
468               DO ji = fs_2, fs_jpim1   ! vector opt.
469                  psi(ji,jj,jk)  = eps(ji,jj,jk)
470               END DO
471            END DO
472         END DO
473         !
474      CASE( 2 )               ! k-w
475         DO jk = 2, jpkm1
476            DO jj = 2, jpjm1
477               DO ji = fs_2, fs_jpim1   ! vector opt.
478                  psi(ji,jj,jk)  = SQRT( en(ji,jj,jk) ) / ( rc0 * mxln(ji,jj,jk) )
479               END DO
480            END DO
481         END DO
482         !
483      CASE( 3 )               ! generic
484         DO jk = 2, jpkm1
485            DO jj = 2, jpjm1
486               DO ji = fs_2, fs_jpim1   ! vector opt.
487                  psi(ji,jj,jk)  = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 
488               END DO
489            END DO
490         END DO
491         !
492      END SELECT
493      !
494      ! Now gls (output in psi)
495      ! -------------------------------
496      ! Resolution of a tridiagonal linear system by a "methode de chasse"
497      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
498      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
499      ! Warning : after this step, en : right hand side of the matrix
500
501      DO jk = 2, jpkm1
502         DO jj = 2, jpjm1
503            DO ji = fs_2, fs_jpim1   ! vector opt.
504               !
505               ! psi / k
506               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
507               !
508               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
509               dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
510               !
511               rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p
512               !
513               ! shear prod. - stratif. destruction
514               prod = rpsi1 * zratio * shear(ji,jj,jk)
515               !
516               ! stratif. destruction
517               buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )
518               !
519               ! shear prod. - stratif. destruction
520               diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
521               !
522               dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
523               !
524               zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term
525               zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
526               !                                                       
527               ! building the matrix
528               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
529               ! lower diagonal
530               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
531                  &                      / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
532               ! upper diagonal
533               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
534                  &                      / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
535               ! diagonal
536               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
537                  &                       + rdt * zdiss * tmask(ji,jj,jk)
538               !
539               ! right hand side in psi
540               psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
541            END DO
542         END DO
543      END DO
544      !
545      z_elem_b(:,:,jpk) = 1._wp
546
547      ! Surface boundary condition on psi
548      ! ---------------------------------
549      !
550      SELECT CASE ( nn_psibc_surf )
551      !
552      CASE ( 0 )             ! Dirichlet boundary conditions
553         !
554         IF( ln_crban ) THEN       ! Wave induced mixing case
555            !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
556            !                      ! balance between the production and the dissipation terms including the wave effect
557            zdep(:,:) = rl_sf * zhsro(:,:)
558            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
559            z_elem_a(:,:,1) = psi(:,:,1)
560            z_elem_c(:,:,1) = 0._wp
561            z_elem_b(:,:,1) = 1._wp
562            !
563            ! one level below
564            zex1 = (rmm*ra_sf+rnn)
565            zex2 = (rmm*ra_sf)
566            zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2
567            psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)
568            z_elem_a(:,:,2) = 0._wp
569            z_elem_c(:,:,2) = 0._wp
570            z_elem_b(:,:,2) = 1._wp
571            !
572         ELSE                   ! No wave induced mixing case
573            !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs
574            !                      ! balance between the production and the dissipation terms
575            !
576            zdep(:,:) = vkarmn * zhsro(:,:)
577            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
578            z_elem_a(:,:,1) = psi(:,:,1)
579            z_elem_c(:,:,1) = 0._wp
580            z_elem_b(:,:,1) = 1._wp
581            !
582            ! one level below
583            zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) )
584            psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
585            z_elem_a(:,:,2) = 0._wp
586            z_elem_c(:,:,2) = 0._wp
587            z_elem_b(:,:,2) = 1.
588            !
589         ENDIF
590         !
591      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
592         !
593         IF( ln_crban ) THEN     ! Wave induced mixing case
594            !
595            zdep(:,:) = rl_sf * zhsro(:,:)
596            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
597            z_elem_a(:,:,1) = psi(:,:,1)
598            z_elem_c(:,:,1) = 0._wp
599            z_elem_b(:,:,1) = 1._wp
600            !
601            ! Neumann condition at k=2
602            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
603            z_elem_a(:,:,2) = 0._wp
604            !
605            ! Set psi vertical flux at the surface:
606            zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)
607            zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
608               &                   * en(:,:,1)**rmm * zdep         
609            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
610            !
611      ELSE                   ! No wave induced mixing
612            !
613            zdep(:,:) = vkarmn * zhsro(:,:)
614            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
615            z_elem_a(:,:,1) = psi(:,:,1)
616            z_elem_c(:,:,1) = 0._wp
617            z_elem_b(:,:,1) = 1._wp
618            !
619            ! Neumann condition at k=2
620            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
621            z_elem_a(ji,jj,2) = 0._wp
622            !
623            ! Set psi vertical flux at the surface:
624            zdep(:,:)  = zhsro(:,:) + fsdept(:,:,1)
625            zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp)
626            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
627            !     
628         ENDIF
629         !
630      END SELECT
631
632      ! Bottom boundary condition on psi
633      ! --------------------------------
634      !
635      SELECT CASE ( nn_psibc_bot )
636      !
637      CASE ( 0 )             ! Dirichlet
638         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro
639         !                      ! Balance between the production and the dissipation terms
640!CDIR NOVERRCHK
641         DO jj = 2, jpjm1
642!CDIR NOVERRCHK
643            DO ji = fs_2, fs_jpim1   ! vector opt.
644               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
645               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
646               zdep(ji,jj) = vkarmn * hbro
647               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
648               z_elem_a(ji,jj,ibot) = 0._wp
649               z_elem_c(ji,jj,ibot) = 0._wp
650               z_elem_b(ji,jj,ibot) = 1._wp
651               !
652               ! Just above last level, Dirichlet condition again (GOTM like)
653               zdep(ji,jj) = vkarmn * ( hbro + fse3t(ji,jj,ibotm1) )
654               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
655               z_elem_a(ji,jj,ibotm1) = 0._wp
656               z_elem_c(ji,jj,ibotm1) = 0._wp
657               z_elem_b(ji,jj,ibotm1) = 1._wp
658            END DO
659         END DO
660         !
661      CASE ( 1 )             ! Neumman boundary condition
662         !                     
663!CDIR NOVERRCHK
664         DO jj = 2, jpjm1
665!CDIR NOVERRCHK
666            DO ji = fs_2, fs_jpim1   ! vector opt.
667               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
668               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
669               !
670               ! Bottom level Dirichlet condition:
671               zdep(ji,jj) = vkarmn * hbro
672               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
673               !
674               z_elem_a(ji,jj,ibot) = 0._wp
675               z_elem_c(ji,jj,ibot) = 0._wp
676               z_elem_b(ji,jj,ibot) = 1._wp
677               !
678               ! Just above last level: Neumann condition with flux injection
679               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b
680               z_elem_c(ji,jj,ibotm1) = 0.
681               !
682               ! Set psi vertical flux at the bottom:
683               zdep(ji,jj) = hbro + 0.5_wp*fse3t(ji,jj,ibotm1)
684               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   &
685                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
686               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)
687            END DO
688         END DO
689         !
690      END SELECT
691
692      ! Matrix inversion
693      ! ----------------
694      !
695      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
696         DO jj = 2, jpjm1
697            DO ji = fs_2, fs_jpim1    ! vector opt.
698               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
699            END DO
700         END DO
701      END DO
702      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
703         DO jj = 2, jpjm1
704            DO ji = fs_2, fs_jpim1    ! vector opt.
705               z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
706            END DO
707         END DO
708      END DO
709      DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
710         DO jj = 2, jpjm1
711            DO ji = fs_2, fs_jpim1    ! vector opt.
712               psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
713            END DO
714         END DO
715      END DO
716
717      ! Set dissipation
718      !----------------
719
720      SELECT CASE ( nn_clos )
721      !
722      CASE( 0 )               ! k-kl  (Mellor-Yamada)
723         DO jk = 1, jpkm1
724            DO jj = 2, jpjm1
725               DO ji = fs_2, fs_jpim1   ! vector opt.
726                  eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk)
727               END DO
728            END DO
729         END DO
730         !
731      CASE( 1 )               ! k-eps
732         DO jk = 1, jpkm1
733            DO jj = 2, jpjm1
734               DO ji = fs_2, fs_jpim1   ! vector opt.
735                  eps(ji,jj,jk) = psi(ji,jj,jk)
736               END DO
737            END DO
738         END DO
739         !
740      CASE( 2 )               ! k-w
741         DO jk = 1, jpkm1
742            DO jj = 2, jpjm1
743               DO ji = fs_2, fs_jpim1   ! vector opt.
744                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
745               END DO
746            END DO
747         END DO
748         !
749      CASE( 3 )               ! generic
750         zcoef = rc0**( 3._wp  + rpp/rnn )
751         zex1  =      ( 1.5_wp + rmm/rnn )
752         zex2  = -1._wp / rnn
753         DO jk = 1, jpkm1
754            DO jj = 2, jpjm1
755               DO ji = fs_2, fs_jpim1   ! vector opt.
756                  eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
757               END DO
758            END DO
759         END DO
760         !
761      END SELECT
762
763      ! Limit dissipation rate under stable stratification
764      ! --------------------------------------------------
765      DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time
766         DO jj = 2, jpjm1
767            DO ji = fs_2, fs_jpim1    ! vector opt.
768               ! limitation
769               eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
770               mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
771               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
772               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
773               mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk)  )
774            END DO
775         END DO
776      END DO 
777
778      !
779      ! Stability function and vertical viscosity and diffusivity
780      ! ---------------------------------------------------------
781      !
782      SELECT CASE ( nn_stab_func )
783      !
784      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
785         DO jk = 2, jpkm1
786            DO jj = 2, jpjm1
787               DO ji = fs_2, fs_jpim1   ! vector opt.
788                  ! zcof =  l²/q²
789                  zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
790                  ! Gh = -N²l²/q²
791                  gh = - rn2(ji,jj,jk) * zcof
792                  gh = MIN( gh, rgh0   )
793                  gh = MAX( gh, rghmin )
794                  ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
795                  sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) )
796                  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)
797                  !
798                  ! Store stability function in avmu and avmv
799                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
800                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
801               END DO
802            END DO
803         END DO
804         !
805      CASE ( 2, 3 )               ! Canuto stability functions
806         DO jk = 2, jpkm1
807            DO jj = 2, jpjm1
808               DO ji = fs_2, fs_jpim1   ! vector opt.
809                  ! zcof =  l²/q²
810                  zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
811                  ! Gh = -N²l²/q²
812                  gh = - rn2(ji,jj,jk) * zcof
813                  gh = MIN( gh, rgh0   )
814                  gh = MAX( gh, rghmin )
815                  gh = gh * rf6
816                  ! Gm =  M²l²/q² Shear number
817                  shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall )
818                  gm = MAX( shr * zcof , 1.e-10 )
819                  gm = gm * rf6
820                  gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm )
821                  ! Stability functions from Canuto
822                  rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm
823                  sm = (rs0 - rs1*gh + rs2*gm) / rcff
824                  sh = (rs4 - rs5*gh + rs6*gm) / rcff
825                  !
826                  ! Store stability function in avmu and avmv
827                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
828                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
829               END DO
830            END DO
831         END DO
832         !
833      END SELECT
834
835      ! Boundary conditions on stability functions for momentum (Neumann):
836      ! Lines below are useless if GOTM style Dirichlet conditions are used
837      zcoef = rcm_sf / SQRT( 2._wp )
838      DO jj = 2, jpjm1
839         DO ji = fs_2, fs_jpim1   ! vector opt.
840            avmv(ji,jj,1) = zcoef
841         END DO
842      END DO
843      zcoef = rc0 / SQRT( 2._wp )
844      DO jj = 2, jpjm1
845         DO ji = fs_2, fs_jpim1   ! vector opt.
846            avmv(ji,jj,mbkt(ji,jj)+1) = zcoef
847         END DO
848      END DO
849
850      ! Compute diffusivities/viscosities
851      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
852      DO jk = 1, jpk
853         DO jj = 2, jpjm1
854            DO ji = fs_2, fs_jpim1   ! vector opt.
855               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)
856               zav           = zsqen * avmu(ji,jj,jk)
857               avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine
858               zav           = zsqen * avmv(ji,jj,jk)
859               avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom
860            END DO
861         END DO
862      END DO
863      !
864      ! Lateral boundary conditions (sign unchanged)
865      avt(:,:,1)  = 0._wp
866      CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
867
868      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
869         DO jj = 2, jpjm1
870            DO ji = fs_2, fs_jpim1   ! vector opt.
871               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
872               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
873            END DO
874         END DO
875      END DO
876      avmu(:,:,1) = 0._wp             ;   avmv(:,:,1) = 0._wp                 ! set surface to zero
877      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )       ! Lateral boundary conditions
878
879      IF(ln_ctl) THEN
880         CALL prt_ctl( tab3d_1=en  , clinfo1=' gls  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
881         CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls  - u: ', mask1=umask,                   &
882            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
883      ENDIF
884      !
885      IF( wrk_not_released(2, 1,2,3)     .OR. &
886          wrk_not_released(3, 1,2,3,4,5)  )   CALL ctl_stop('zdf_gls: failed to release workspace arrays')
887      !
888   END SUBROUTINE zdf_gls
889
890
891   SUBROUTINE zdf_gls_init
892      !!----------------------------------------------------------------------
893      !!                  ***  ROUTINE zdf_gls_init  ***
894      !!                     
895      !! ** Purpose :   Initialization of the vertical eddy diffivity and
896      !!      viscosity when using a gls turbulent closure scheme
897      !!
898      !! ** Method  :   Read the namzdf_gls namelist and check the parameters
899      !!      called at the first timestep (nit000)
900      !!
901      !! ** input   :   Namlist namzdf_gls
902      !!
903      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
904      !!
905      !!----------------------------------------------------------------------
906      USE dynzdf_exp
907      USE trazdf_exp
908      !
909      INTEGER ::   jk    ! dummy loop indices
910      REAL(wp)::   zcr   ! local scalar
911      !!
912      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
913         &            rn_clim_galp, ln_crban, ln_sigpsi,     &
914         &            rn_crban, rn_charn,                    &
915         &            nn_tkebc_surf, nn_tkebc_bot,           &
916         &            nn_psibc_surf, nn_psibc_bot,           &
917         &            nn_stab_func, nn_clos
918      !!----------------------------------------------------------
919
920      REWIND( numnam )                 !* Read Namelist namzdf_gls
921      READ  ( numnam, namzdf_gls )
922
923      IF(lwp) THEN                     !* Control print
924         WRITE(numout,*)
925         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
926         WRITE(numout,*) '~~~~~~~~~~~~'
927         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
928         WRITE(numout,*) '      minimum value of en                           rn_emin       = ', rn_emin
929         WRITE(numout,*) '      minimum value of eps                          rn_epsmin     = ', rn_epsmin
930         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim = ', ln_length_lim
931         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp  = ', rn_clim_galp
932         WRITE(numout,*) '      TKE Surface boundary condition                nn_tkebc_surf = ', nn_tkebc_surf
933         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_tkebc_bot  = ', nn_tkebc_bot
934         WRITE(numout,*) '      PSI Surface boundary condition                nn_psibc_surf = ', nn_psibc_surf
935         WRITE(numout,*) '      PSI Bottom boundary condition                 nn_psibc_bot  = ', nn_psibc_bot
936         WRITE(numout,*) '      Craig and Banner scheme                       ln_crban      = ', ln_crban
937         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi     = ', ln_sigpsi
938         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
939         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
940         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
941         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
942         WRITE(numout,*) '   Hard coded parameters'
943         WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro
944         WRITE(numout,*) '      Bottom roughness (m)                          hbro          = ', hbro
945      ENDIF
946
947      !                                !* allocate gls arrays
948      IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
949
950      !                                !* Check of some namelist values
951      IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' )
952      IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' )
953      IF( nn_tkebc_bot  < 0 .OR. nn_tkebc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' )
954      IF( nn_psibc_bot  < 0 .OR. nn_psibc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' )
955      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
956      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
957
958      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
959      !
960      CASE( 0 )                              ! k-kl  (Mellor-Yamada)
961         !
962         IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
963         rpp     = 0._wp
964         rmm     = 1._wp
965         rnn     = 1._wp
966         rsc_tke = 1.96_wp
967         rsc_psi = 1.96_wp
968         rpsi1   = 0.9_wp
969         rpsi3p  = 1._wp
970         rpsi2   = 0.5_wp
971         !
972         SELECT CASE ( nn_stab_func )
973         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
974         CASE( 2 )      ;   rpsi3m = 2.38_wp       ! Canuto A stability functions
975         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
976         END SELECT
977         !
978      CASE( 1 )                              ! k-eps
979         !
980         IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
981         rpp     =  3._wp
982         rmm     =  1.5_wp
983         rnn     = -1._wp
984         rsc_tke =  1._wp
985         rsc_psi =  1.3_wp  ! Schmidt number for psi
986         rpsi1   =  1.44_wp
987         rpsi3p  =  1._wp
988         rpsi2   =  1.92_wp
989         !
990         SELECT CASE ( nn_stab_func )
991         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
992         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
993         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
994         END SELECT
995         !
996      CASE( 2 )                              ! k-omega
997         !
998         IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
999         rpp     = -1._wp
1000         rmm     =  0.5_wp
1001         rnn     = -1._wp
1002         rsc_tke =  2._wp
1003         rsc_psi =  2._wp
1004         rpsi1   =  0.555_wp
1005         rpsi3p  =  1._wp
1006         rpsi2   =  0.833_wp
1007         !
1008         SELECT CASE ( nn_stab_func )
1009         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
1010         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
1011         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
1012         END SELECT
1013         !
1014      CASE( 3 )                              ! generic
1015         !
1016         IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
1017         rpp     = 2._wp
1018         rmm     = 1._wp
1019         rnn     = -0.67_wp
1020         rsc_tke = 0.8_wp
1021         rsc_psi = 1.07_wp
1022         rpsi1   = 1._wp
1023         rpsi3p  = 1._wp
1024         rpsi2   = 1.22_wp
1025         !
1026         SELECT CASE ( nn_stab_func )
1027         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
1028         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
1029         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
1030         END SELECT
1031         !
1032      END SELECT
1033
1034      !
1035      SELECT CASE ( nn_stab_func )     !* set the parameters of the stability functions
1036      !
1037      CASE ( 0 )                             ! Galperin stability functions
1038         !
1039         IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
1040         rc2     =  0._wp
1041         rc3     =  0._wp
1042         rc_diff =  1._wp
1043         rc0     =  0.5544_wp
1044         rcm_sf  =  0.9884_wp
1045         rghmin  = -0.28_wp
1046         rgh0    =  0.0233_wp
1047         rghcri  =  0.02_wp
1048         !
1049      CASE ( 1 )                             ! Kantha-Clayson stability functions
1050         !
1051         IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
1052         rc2     =  0.7_wp
1053         rc3     =  0.2_wp
1054         rc_diff =  1._wp
1055         rc0     =  0.5544_wp
1056         rcm_sf  =  0.9884_wp
1057         rghmin  = -0.28_wp
1058         rgh0    =  0.0233_wp
1059         rghcri  =  0.02_wp
1060         !
1061      CASE ( 2 )                             ! Canuto A stability functions
1062         !
1063         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
1064         rs0 = 1.5_wp * rl1 * rl5*rl5
1065         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
1066         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
1067         rs4 = 2._wp * rl5
1068         rs5 = 2._wp * rl4
1069         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
1070            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
1071         rd0 = 3._wp * rl5*rl5
1072         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
1073         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
1074         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
1075         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
1076         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
1077         rc0 = 0.5268_wp
1078         rf6 = 8._wp / (rc0**6._wp)
1079         rc_diff = SQRT(2._wp) / (rc0**3._wp)
1080         rcm_sf  =  0.7310_wp
1081         rghmin  = -0.28_wp
1082         rgh0    =  0.0329_wp
1083         rghcri  =  0.03_wp
1084         !
1085      CASE ( 3 )                             ! Canuto B stability functions
1086         !
1087         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
1088         rs0 = 1.5_wp * rm1 * rm5*rm5
1089         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
1090         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
1091         rs4 = 2._wp * rm5
1092         rs5 = 2._wp * rm4
1093         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)
1094         rd0 = 3._wp * rm5*rm5
1095         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
1096         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
1097         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
1098         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
1099         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
1100         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
1101         rf6 = 8._wp / ( rc0**6._wp )
1102         rc_diff = SQRT(2._wp)/(rc0**3.)
1103         rcm_sf  =  0.7470_wp
1104         rghmin  = -0.28_wp
1105         rgh0    =  0.0444_wp
1106         rghcri  =  0.0414_wp
1107         !
1108      END SELECT
1109   
1110      !                                !* Set Schmidt number for psi diffusion in the wave breaking case
1111      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1112      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1113      IF( ln_sigpsi .AND. ln_crban ) THEN
1114         zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn
1115         rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf )                       & 
1116        &         * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn   &
1117        &           + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm )                                 
1118      ELSE
1119         rsc_psi0 = rsc_psi
1120      ENDIF
1121 
1122      !                                !* Shear free turbulence parameters
1123      !
1124      ra_sf  = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke )   &
1125         &                                      - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
1126      rl_sf  = rc0 * SQRT( rc0 / rcm_sf )                                                                   &
1127         &         * SQRT(  (  (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke                                &
1128         &                   + 12._wp * rsc_psi0 * rpsi2                                                    &
1129         &                   - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) )  )    &
1130         &                / ( 12._wp*rnn*rnn )                                                              )
1131
1132      !
1133      IF(lwp) THEN                     !* Control print
1134         WRITE(numout,*)
1135         WRITE(numout,*) 'Limit values'
1136         WRITE(numout,*) '~~~~~~~~~~~~'
1137         WRITE(numout,*) 'Parameter  m = ',rmm
1138         WRITE(numout,*) 'Parameter  n = ',rnn
1139         WRITE(numout,*) 'Parameter  p = ',rpp
1140         WRITE(numout,*) 'rpsi1   = ',rpsi1
1141         WRITE(numout,*) 'rpsi2   = ',rpsi2
1142         WRITE(numout,*) 'rpsi3m  = ',rpsi3m
1143         WRITE(numout,*) 'rpsi3p  = ',rpsi3p
1144         WRITE(numout,*) 'rsc_tke = ',rsc_tke
1145         WRITE(numout,*) 'rsc_psi = ',rsc_psi
1146         WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0
1147         WRITE(numout,*) 'rc0     = ',rc0
1148         WRITE(numout,*)
1149         WRITE(numout,*) 'Shear free turbulence parameters:'
1150         WRITE(numout,*) 'rcm_sf  = ',rcm_sf
1151         WRITE(numout,*) 'ra_sf   = ',ra_sf
1152         WRITE(numout,*) 'rl_sf   = ',rl_sf
1153         WRITE(numout,*)
1154      ENDIF
1155
1156      !                                !* Constants initialization
1157      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
1158      rc03  = rc02 * rc0
1159      rc04  = rc03 * rc0
1160      rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp
1161      rsbc_mb   = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp)               ! Surf. bound. cond. from Mellor and Blumberg
1162      rsbc_std  = 3.75_wp                                                  ! Surf. bound. cond. standard (prod=diss)
1163      rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.  Dirichlet + Wave breaking
1164      rsbc_tke2 = 0.5_wp / rau0
1165      rsbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking
1166      rsbc_zs   = rn_charn / grav                                                        ! Charnock formula
1167      rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn                           ! Dirichlet + Wave breaking
1168      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking
1169      rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + Wave breaking
1170      rfact_tke = -0.5_wp / rsc_tke * rdt               ! Cst used for the Diffusion term of tke
1171      rfact_psi = -0.5_wp / rsc_psi * rdt               ! Cst used for the Diffusion term of tke
1172
1173      !                                !* Wall proximity function
1174      zwall (:,:,:) = 1._wp * tmask(:,:,:)
1175
1176      !                                !* set vertical eddy coef. to the background value
1177      DO jk = 1, jpk
1178         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1179         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
1180         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1181         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1182      END DO
1183      !                             
1184      CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
1185      !
1186   END SUBROUTINE zdf_gls_init
1187
1188
1189   SUBROUTINE gls_rst( kt, cdrw )
1190      !!---------------------------------------------------------------------
1191      !!                   ***  ROUTINE ts_rst  ***
1192      !!                     
1193      !! ** Purpose :   Read or write TKE file (en) in restart file
1194      !!
1195      !! ** Method  :   use of IOM library
1196      !!                if the restart does not contain TKE, en is either
1197      !!                set to rn_emin or recomputed (nn_igls/=0)
1198      !!----------------------------------------------------------------------
1199      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1200      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1201      !
1202      INTEGER ::   jit, jk   ! dummy loop indices
1203      INTEGER ::   id1, id2, id3, id4, id5, id6, id7, id8
1204      INTEGER ::   ji, jj, ikbu, ikbv
1205      REAL(wp)::   cbx, cby
1206      !!----------------------------------------------------------------------
1207      !
1208      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1209         !                                   ! ---------------
1210         IF( ln_rstart ) THEN                   !* Read the restart file
1211            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
1212            id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
1213            id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
1214            id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
1215            id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
1216            id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
1217            !
1218            IF( MIN( id1, id2, id3, id4, id5, id6, id7, id8 ) > 0 ) THEN        ! all required arrays exist
1219               CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
1220               CALL iom_get( numror, jpdom_autoglo, 'avt'   , avt    )
1221               CALL iom_get( numror, jpdom_autoglo, 'avm'   , avm    )
1222               CALL iom_get( numror, jpdom_autoglo, 'avmu'  , avmu   )
1223               CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   )
1224               CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
1225            ELSE                       
1226               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
1227               en  (:,:,:) = rn_emin
1228               mxln(:,:,:) = 0.001       
1229               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
1230            ENDIF
1231         ELSE                                   !* Start from rest
1232            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
1233            en  (:,:,:) = rn_emin
1234            mxln(:,:,:) = 0.001       
1235         ENDIF
1236         !
1237      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1238         !                                   ! -------------------
1239         IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1240         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
1241         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt   )
1242         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm   )
1243         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu  )
1244         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  )
1245         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln  )
1246         !
1247      ENDIF
1248      !
1249   END SUBROUTINE gls_rst
1250
1251#else
1252   !!----------------------------------------------------------------------
1253   !!   Dummy module :                                        NO TKE scheme
1254   !!----------------------------------------------------------------------
1255   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfgls = .FALSE.   !: TKE flag
1256CONTAINS
1257   SUBROUTINE zdf_gls_init           ! Empty routine
1258      WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
1259   END SUBROUTINE zdf_gls_init
1260   SUBROUTINE zdf_gls( kt )          ! Empty routine
1261      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
1262   END SUBROUTINE zdf_gls
1263   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
1264      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1265      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1266      WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw
1267   END SUBROUTINE gls_rst
1268#endif
1269
1270   !!======================================================================
1271END MODULE zdfgls
1272
Note: See TracBrowser for help on using the repository browser.