source: branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9353

Last change on this file since 9353 was 9353, checked in by jchanut, 2 years ago

Tilde coordinate update: first draft

  • Property svn:keywords set to Id
File size: 148.5 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
9   !!                                          vvl option includes z_star and z_tilde coordinates
10   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
11   !!                 !  2018-01  (J. Chanut) improve ztilde robustness
12   !!----------------------------------------------------------------------
13   !!   'key_vvl'                              variable volume
14   !!----------------------------------------------------------------------
15   !!----------------------------------------------------------------------
16   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
17   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
18   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
19   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
20   !!   dom_vvl_rst      : read/write restart file
21   !!   dom_vvl_ctl      : Check the vvl options
22   !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors
23   !!                    : to account for manual changes to e[1,2][u,v] in some Straits
24   !!----------------------------------------------------------------------
25   !! * Modules used
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sbc_oce         ! ocean surface boundary condition
29   USE in_out_manager  ! I/O manager
30   USE iom             ! I/O manager library
31   USE restart         ! ocean restart
32   USE lib_mpp         ! distributed memory computing library
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE wrk_nemo        ! Memory allocation
35   USE timing          ! Timing
36   USE bdy_oce         ! ocean open boundary conditions
37
38   IMPLICIT NONE
39   PRIVATE
40
41   !! * Routine accessibility
42   PUBLIC  dom_vvl_init       ! called by domain.F90
43   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
44   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
45   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
46   PRIVATE dom_vvl_orca_fix   ! called by dom_vvl_interpol
47
48   !!* Namelist nam_vvl
49   LOGICAL , PUBLIC                                      :: ln_vvl_zstar = .FALSE.              ! zstar  vertical coordinate
50   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.             ! ztilde vertical coordinate
51   LOGICAL , PUBLIC                                      :: ln_vvl_layer = .FALSE.              ! level  vertical coordinate
52   LOGICAL                                               :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
53   LOGICAL                                               :: ln_vvl_zstar_at_eqtor = .FALSE.     ! revert to zstar at equator
54   LOGICAL                                               :: ln_vvl_zstar_on_shelf =.FALSE.      ! revert to zstar on shelves
55   LOGICAL                                               :: ln_vvl_kepe =.FALSE.                ! kinetic/potential energy transfer
56   LOGICAL                                               :: ln_vvl_adv_fct =.FALSE.             ! Centred thickness advection
57   LOGICAL                                               :: ln_vvl_adv_cn2 =.TRUE.              ! FCT thickness advection
58   LOGICAL                                               :: ln_vvl_dbg = .FALSE.                ! debug control prints
59   LOGICAL                                               :: ln_vvl_lap                          ! Laplacian thickness diffusion
60   LOGICAL                                               :: ln_vvl_blp                          ! Bilaplacian thickness diffusion
61   LOGICAL                                               :: ln_vvl_regrid                       ! ensure layer separation
62   !                                                                                            ! conservation: not used yet
63   REAL(wp)                                              :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian)
64   REAL(wp)                                              :: rn_ahe3_blp               ! thickness diffusion coefficient (Bilaplacian)
65   REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
66   REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
67   REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
68
69   !! * Module variables
70   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport
71   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn_lf,  un_lf,  vn_lf  ! 1st order filtered arrays         
72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n  ! baroclinic scale factors
73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors
74   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t               ! restoring period for scale factors
75   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv               ! restoring period for low freq. divergence
76   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: tildemask                 ! mask tilde tendency
77   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: hsm, dsm                  !
78   INTEGER         , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: i_int_bot
79
80   !! * Substitutions
81#  include "domzgr_substitute.h90"
82#  include "vectopt_loop_substitute.h90"
83   !!----------------------------------------------------------------------
84   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
85   !! $Id$
86   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
87   !!----------------------------------------------------------------------
88
89CONTAINS
90
91   INTEGER FUNCTION dom_vvl_alloc()
92      !!----------------------------------------------------------------------
93      !!                ***  FUNCTION dom_vvl_alloc  ***
94      !!----------------------------------------------------------------------
95      IF( ln_vvl_zstar ) dom_vvl_alloc = 0
96      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
97         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
98            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
99            &      i_int_bot(jpi,jpj), tildemask(jpi,jpj), hsm(jpi,jpj), dsm(jpi,jpj), &
100            &      STAT = dom_vvl_alloc        )
101         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
102         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
103         un_td = 0.0_wp
104         vn_td = 0.0_wp
105      ENDIF
106      IF( ln_vvl_ztilde ) THEN
107         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj),  &
108            &        hdivn_lf(jpi,jpj,jpk),  &
109            &           un_lf(jpi,jpj,jpk),  &
110            &           vn_lf(jpi,jpj,jpk),  STAT= dom_vvl_alloc )
111         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
112         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
113      ENDIF
114
115   END FUNCTION dom_vvl_alloc
116
117
118   SUBROUTINE dom_vvl_init
119      !!----------------------------------------------------------------------
120      !!                ***  ROUTINE dom_vvl_init  ***
121      !!                   
122      !! ** Purpose :  Initialization of all scale factors, depths
123      !!               and water column heights
124      !!
125      !! ** Method  :  - use restart file and/or initialize
126      !!               - interpolate scale factors
127      !!
128      !! ** Action  : - fse3t_(n/b) and tilde_e3t_(n/b)
129      !!              - Regrid: fse3(u/v)_n
130      !!                        fse3(u/v)_b       
131      !!                        fse3w_n           
132      !!                        fse3(u/v)w_b     
133      !!                        fse3(u/v)w_n     
134      !!                        fsdept_n, fsdepw_n and fsde3w_n
135      !!              - h(t/u/v)_0
136      !!              - frq_rst_e3t and frq_rst_hdv
137      !!
138      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
139      !!----------------------------------------------------------------------
140      USE phycst,  ONLY : rpi, rsmall, rad
141      !! * Local declarations
142      INTEGER ::   ji, jj, jk
143      INTEGER ::   ii0, ii1, ij0, ij1
144      REAL(wp)::   zcoef, zwgt, ztmp, zhmin, zhmax
145      !!----------------------------------------------------------------------
146      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
147
148      IF(lwp) WRITE(numout,*)
149      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
150      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
151
152      ! choose vertical coordinate (z_star, z_tilde or layer)
153      ! ==========================
154      CALL dom_vvl_ctl
155
156      ! Allocate module arrays
157      ! ======================
158      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
159
160      ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
161      ! ============================================================================
162      CALL dom_vvl_rst( nit000, 'READ' )
163      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
164
165      ! Reconstruction of all vertical scale factors at now and before time steps
166      ! =============================================================================
167      ! Horizontal scale factor interpolations
168      ! --------------------------------------
169      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
170      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
171      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
172      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
173      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' )
174      ! Vertical scale factor interpolations
175      ! ------------------------------------
176      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
177      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
178      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
179      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
180      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
181      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
182      ! t- and w- points depth
183      ! ----------------------
184      ! set the isf depth as it is in the initial step
185      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
186      fsdepw_n(:,:,1) = 0.0_wp
187      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
188      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
189      fsdepw_b(:,:,1) = 0.0_wp
190
191      DO jk = 2, jpk
192         DO jj = 1,jpj
193            DO ji = 1,jpi
194              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
195                                                     ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
196                                                     ! 0.5 where jk = mikt 
197               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
198               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
199               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
200                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
201               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
202               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
203               fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  &
204                   &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk)) 
205            END DO
206         END DO
207      END DO
208
209      ! Before depth and Inverse of the local depth of the water column at u- and v- points
210      ! -----------------------------------------------------------------------------------
211      hu_b(:,:) = 0.
212      hv_b(:,:) = 0.
213      DO jk = 1, jpkm1
214         hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
215         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
216      END DO
217      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
218      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
219
220      ! Restoring frequencies for z_tilde coordinate
221      ! ============================================
222      tildemask(:,:) = 1._wp
223
224      IF( ln_vvl_ztilde ) THEN
225         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
226         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
227         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
228         ! tendency mask:
229         !
230         IF( ln_vvl_ztilde_as_zstar ) THEN
231            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
232            frq_rst_e3t(:,:) = 0.0_wp 
233            frq_rst_hdv(:,:) = 1.0_wp / rdt
234            tildemask(:,:) = 0._wp
235         ENDIF
236         
237         IF ( ln_vvl_zstar_at_eqtor ) THEN
238            DO jj = 1, jpj
239               DO ji = 1, jpi
240                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
241                     ! values outside the equatorial band and transition zone (ztilde)
242                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
243!                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
244             
245                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
246                     ! values inside the equatorial band (ztilde as zstar)
247                     frq_rst_e3t(ji,jj) =  0.0_wp
248!                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
249                     tildemask(ji,jj) = 0._wp
250                  ELSE
251                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
252                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
253                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
254                        &                                          * 180._wp / 3.5_wp ) )
255!                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
256!                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
257!                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
258!                        &                                          * 180._wp / 3.5_wp ) )
259                     tildemask(ji,jj) = 0.5_wp * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
260                        &                                                 * 180._wp / 3.5_wp ) )
261                  ENDIF
262               END DO
263            END DO
264            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
265               ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
266               ij0 = 128   ;   ij1 = 135   ;   
267               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
268               tildemask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )   = 0.0_wp
269!               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
270            ENDIF
271         ENDIF
272         !
273         IF ( ln_vvl_zstar_on_shelf ) THEN
274            zhmin =  50._wp
275            zhmax = 100._wp
276            DO jj = 1, jpj
277               DO ji = 1, jpi
278                  zwgt = 1._wp
279                  IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
280                     zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
281                  ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
282                     zwgt = 0._wp
283                  ENDIF
284                  frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt)
285                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt)
286               END DO
287            END DO
288         ENDIF
289         !
290         ztmp = MAXVAL( frq_rst_hdv(:,:) )
291         IF( lk_mpp )   CALL mpp_max( ztmp )                 ! max over the global domain
292         !
293         IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' )
294         !
295      ENDIF
296
297      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
298
299   END SUBROUTINE dom_vvl_init
300
301
302   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
303      !!----------------------------------------------------------------------
304      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
305      !!                   
306      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
307      !!                 tranxt and dynspg routines
308      !!
309      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
310      !!               - z_tilde_case: after scale factor increment =
311      !!                                    high frequency part of horizontal divergence
312      !!                                  + retsoring towards the background grid
313      !!                                  + thickness difusion
314      !!                               Then repartition of ssh INCREMENT proportionnaly
315      !!                               to the "baroclinic" level thickness.
316      !!
317      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
318      !!               - tilde_e3t_a: after increment of vertical scale factor
319      !!                              in z_tilde case
320      !!               - fse3(t/u/v)_a
321      !!
322      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
323      !!----------------------------------------------------------------------
324      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t, ztu, ztv
325      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
326      !! * Arguments
327      INTEGER, INTENT( in )                  :: kt                    ! time step
328      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence
329      !! * Local declarations
330      INTEGER                                :: ji, jj, jk            ! dummy loop indices
331      INTEGER                                :: ib, ib_bdy, ip, jp    !   "     "     "
332      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
333      REAL(wp)                               :: z2dt                  ! temporary scalars
334      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
335      REAL(wp)                               :: zalpha, zwgt          ! temporary scalars
336      REAL(wp)                               :: zdu, zdv
337      LOGICAL                                :: ll_do_bclinic         ! temporary logical
338      !!----------------------------------------------------------------------
339      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
340      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
341      CALL wrk_alloc( jpi, jpj, jpk, ze3t, ztu, ztv )
342
343      IF(kt == nit000)   THEN
344         IF(lwp) WRITE(numout,*)
345         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
346         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
347      ENDIF
348
349      ll_do_bclinic = .TRUE.
350      IF( PRESENT(kcall) ) THEN
351         IF ( kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE.
352      ENDIF
353
354      ! ******************************* !
355      ! After acale factors at t-points !
356      ! ******************************* !
357
358      !                                                ! --------------------------------------------- !
359                                                       ! z_star coordinate and barotropic z-tilde part !
360      !                                                ! --------------------------------------------- !
361
362      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - ssmask(:,:) )
363      DO jk = 1, jpkm1
364         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
365         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
366      END DO
367     
368      IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
369         
370         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
371         un_td(:,:,:) = 0.0_wp        ! Transport corrections
372         vn_td(:,:,:) = 0.0_wp
373
374         zhdiv(:,:) = 0.
375         DO jk = 1, jpkm1
376            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
377         END DO
378         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )     
379
380         ! Thickness advection:
381         ! --------------------
382         ! Set advection velocities and source term
383         IF ( ln_vvl_ztilde ) THEN
384            !
385            IF ((kt==nit000).AND.(neuler==0)) THEN
386               DO jk = 1, jpkm1
387                  ztu(:,:,jk) = un(:,:,jk)
388                  ztv(:,:,jk) = vn(:,:,jk)
389               END DO
390            ELSE
391               DO jk = 1, jpkm1
392                  ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk)
393                  ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk)
394                  tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk)
395               END DO
396            ENDIF
397            !
398         ELSEIF ( ln_vvl_layer ) THEN
399            !
400            DO jk = 1, jpkm1
401               ztu(:,:,jk) = un(:,:,jk)
402               ztv(:,:,jk) = vn(:,:,jk)
403            END DO
404            !
405         ENDIF
406         !
407         ! Do advection
408         IF     (ln_vvl_adv_fct) THEN
409            CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv )
410            !
411         ELSEIF (ln_vvl_adv_cn2) THEN
412            DO jk = 1, jpkm1
413               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * hdivn(:,:,jk)
414            END DO
415         ENDIF
416         !
417         ! Thickness anomlaly diffusion:
418         ! -----------------------------
419         zwu(:,:) = 0.0_wp
420         zwv(:,:) = 0.0_wp
421         ztu(:,:,:) = 0.0_wp
422         ztv(:,:,:) = 0.0_wp
423
424         IF ( ln_vvl_blp ) THEN  ! Bilaplacian
425            DO jk = 1, jpkm1
426               DO jj = 1, jpjm1                 ! First derivative (gradient)
427                  DO ji = 1, fs_jpim1   ! vector opt.                 
428                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
429                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
430                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
431                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
432                  END DO
433               END DO
434
435               DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
436                  DO ji = fs_2, fs_jpim1 ! vector opt.
437                     zht(ji,jj) =  rn_ahe3_blp * r1_e12t(ji,jj) * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
438                        &                                           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
439                  END DO
440               END DO
441
442#if defined key_bdy
443               DO ib_bdy=1, nb_bdy
444                  DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
445                     ji = idx_bdy(ib_bdy)%nbi(ib,1)
446                     jj = idx_bdy(ib_bdy)%nbj(ib,1)
447                     zht(ji,jj) = 0._wp
448                  END DO
449               END DO
450#endif
451               CALL lbc_lnk( zht, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
452
453               DO jj = 1, jpjm1                 ! third derivative (gradient)
454                  DO ji = 1, fs_jpim1   ! vector opt.
455                     ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) )
456                     ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( zht(ji  ,jj+1) - zht(ji,jj) )
457                     zwu(ji,jj) = zwu(ji,jj) + ztu(ji,jj,jk)
458                     zwv(ji,jj) = zwv(ji,jj) + ztv(ji,jj,jk)
459                  END DO
460               END DO
461            END DO
462         ENDIF
463
464         IF ( ln_vvl_lap ) THEN    ! Laplacian
465            DO jk = 1, jpkm1                    ! First derivative (gradient)
466               DO jj = 1, jpjm1
467                  DO ji = 1, fs_jpim1   ! vector opt.                 
468                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
469                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
470                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
471                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
472                     zwu(ji,jj) = zwu(ji,jj) + zdu
473                     zwv(ji,jj) = zwv(ji,jj) + zdv
474                     ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu
475                     ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv
476                  END DO
477               END DO
478            END DO
479         ENDIF 
480
481         ! Ensure barotropic fluxes are null:
482!         DO jj = 1, jpj
483!            DO ji = 1, jpi
484!               DO jk = 1, jpkm1
485!                  ztu(ji,jj,jk) = ztu(ji,jj,jk) - zwu(ji,jj)*fse3u_b(ji,jj,jk)*hur_b(ji,jj)*umask(ji,jj,jk)
486!                  ztv(ji,jj,jk) = ztv(ji,jj,jk) - zwv(ji,jj)*fse3v_b(ji,jj,jk)*hvr_b(ji,jj)*vmask(ji,jj,jk)
487!               END DO
488!            END DO
489!         END DO
490         DO jj = 1, jpj
491            DO ji = 1, jpi
492               ztu(ji,jj,mbku(ji,jj)) = ztu(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
493               ztv(ji,jj,mbkv(ji,jj)) = ztv(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
494            END DO
495         END DO
496
497         ! divergence of diffusive fluxes
498         DO jk = 1, jpkm1
499            DO jj = 2, jpjm1
500               DO ji = fs_2, fs_jpim1   ! vector opt.
501                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk) - ztu(ji,jj,jk)    &
502                     &                                          +     ztv(ji  ,jj-1,jk) - ztv(ji,jj,jk)    &
503                     &                                            ) * r1_e12t(ji,jj)
504               END DO
505            END DO
506         END DO
507 
508         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:)
509         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:)
510         CALL lbc_lnk( un_td , 'U' , -1.)
511         CALL lbc_lnk( vn_td , 'V' , -1.) 
512         !
513         !
514         ! Restoring:
515         IF( ln_vvl_ztilde ) THEN
516            DO jk = 1, jpk
517               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
518            END DO
519         ENDIF
520
521         ! Remove "external thickness" tendency:
522         DO jk = 1, jpkm1
523            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   fse3t_n(:,:,jk) * zhdiv(:,:)
524         END DO 
525                   
526         ! Leapfrog time stepping
527         ! ~~~~~~~~~~~~~~~~~~~~~~
528         IF( neuler == 0 .AND. kt == nit000 ) THEN
529            z2dt =  rdt
530         ELSE
531            z2dt = 2.0_wp * rdt
532         ENDIF
533
534         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
535
536         ! Revert to zstar locally:
537         ! ~~~~~~~~~~~~~~~~~~~~~~~~
538         DO jk=1,jpkm1
539            tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) * tildemask(:,:) 
540         END DO
541
542         ! Ensure layer separation:
543         ! ~~~~~~~~~~~~~~~~~~~~~~~~
544         IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt )
545 
546         ! Boundary conditions:
547         ! ~~~~~~~~~~~~~~~~~~~~
548#if defined key_bdy
549         DO ib_bdy=1, nb_bdy
550            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
551!!            DO ib = 1, idx_bdy(ib_bdy)%nblen(1)
552               ji = idx_bdy(ib_bdy)%nbi(ib,1)
553               jj = idx_bdy(ib_bdy)%nbj(ib,1)
554               zwgt = idx_bdy(ib_bdy)%nbw(ib,1)
555               ip = bdytmask(ji+1,jj  ) - bdytmask(ji-1,jj  )
556               jp = bdytmask(ji  ,jj+1) - bdytmask(ji  ,jj-1)
557               DO jk = 1, jpkm1 
558                  tilde_e3t_a(ji,jj,jk) = 0.e0
559!!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt)
560!!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk)
561               END DO
562            END DO
563         END DO
564#endif
565         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )       
566
567         ! Maximum deformation control
568         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
569         ze3t(:,:,jpk) = 0.0_wp
570         DO jk = 1, jpkm1
571            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
572         END DO
573         z_tmax = MAXVAL( ze3t(:,:,:) )
574         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
575         z_tmin = MINVAL( ze3t(:,:,:) )
576         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
577         ! - ML - test: for the moment, stop simulation for too large e3_t variations
578         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
579            IF( lk_mpp ) THEN
580               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
581               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
582            ELSE
583               ijk_max = MAXLOC( ze3t(:,:,:) )
584               ijk_max(1) = ijk_max(1) + nimpp - 1
585               ijk_max(2) = ijk_max(2) + njmpp - 1
586               ijk_min = MINLOC( ze3t(:,:,:) )
587               ijk_min(1) = ijk_min(1) + nimpp - 1
588               ijk_min(2) = ijk_min(2) + njmpp - 1
589            ENDIF
590            IF (lwp) THEN
591               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
592               WRITE(numout, *) 'at i, j, k=', ijk_max
593               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
594               WRITE(numout, *) 'at i, j, k=', ijk_min           
595               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
596            ENDIF
597         ENDIF
598
599         IF ( ln_vvl_ztilde ) THEN
600            zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
601            DO jk = 1, jpkm1
602               ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:) + un_td(:,:,jk)
603               ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk)
604               ze3t(:,:,jk) = -fse3t_n(:,:,jk) * zhdiv(:,:)
605!
606!               DO jj = 2, jpjm1
607!                  DO ji = fs_2, fs_jpim1   ! vector opt.
608!                     
609!                     ze3t(ji,jj,jk) =   -fse3t_n(ji,jj,jk) * zhdiv(ji,jj)   &
610!                        & + (  un_td(ji,jj,jk) - un_td(ji-1,jj  ,jk)        &
611!                        &    + vn_td(ji,jj,jk) - vn_td(ji  ,jj-1,jk)  )     &
612!                        & / ( e1t(ji,jj) * e2t(ji,jj) )
613!                  END DO
614!               END DO
615            END DO
616            !
617            un_lf(:,:,:)    = un_lf(:,:,:)    * (1._wp - zalpha) + zalpha * ztu(:,:,:)
618            vn_lf(:,:,:)    = vn_lf(:,:,:)    * (1._wp - zalpha) + zalpha * ztv(:,:,:)
619            hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
620         ENDIF
621
622      ENDIF
623
624      IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.(.NOT.ll_do_bclinic) ) THEN
625         zhdiv(:,:) = 0.
626         DO jk = 1, jpkm1
627            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * (hdivn(:,:,jk) - hdivb(:,:,jk))
628         END DO
629         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
630
631         IF( neuler == 0 .AND. kt == nit000 ) THEN
632            z2dt =  rdt
633         ELSE
634            z2dt = 2.0_wp * rdt
635         ENDIF
636
637         DO jk = 1, jpkm1
638            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * & 
639                                  & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:))
640         END DO
641      ENDIF
642
643      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
644      !                                             ! ---baroclinic part--------- !
645         DO jk = 1, jpkm1
646            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                     
647         END DO
648      ENDIF
649
650      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
651         !
652         zht(:,:) = 0.0_wp
653         DO jk = 1, jpkm1
654            zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
655         END DO
656         IF( lwp ) WRITE(numout, *) 'kt =', kt
657         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
658            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 )
659            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
660            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
661         END IF
662         !
663         z_tmin = MINVAL( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0  )
664         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
665         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin         
666         !
667         z_tmin = MINVAL( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 )
668         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
669         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3u_n) =', z_tmin
670         !
671         z_tmin = MINVAL( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 )
672         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
673         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3v_n) =', z_tmin
674         !
675         z_tmin = MINVAL( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 )
676         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
677         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3f_n) =', z_tmin
678         !
679         zht(:,:) = 0.0_wp
680         DO jk = 1, jpkm1
681            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
682         END DO
683         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
684         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
685         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn  -SUM(fse3t_n))) =', z_tmax
686         !
687         zht(:,:) = 0.0_wp
688         DO jk = 1, jpkm1
689            zht(:,:) = zht(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
690         END DO
691         zwu(:,:) = 0._wp
692         DO jj = 1, jpjm1
693            DO ji = 1, fs_jpim1   ! vector opt.
694               zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e12u(ji,jj)                             &
695                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji+1,jj) * sshn(ji+1,jj) )
696            END DO
697         END DO
698         CALL lbc_lnk( zwu(:,:), 'U', 1._wp )
699         z_tmax = MAXVAL( umask(:,:,1) * umask_i(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) )
700         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
701         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(fse3u_n))) =', z_tmax
702         !
703         zht(:,:) = 0.0_wp
704         DO jk = 1, jpkm1
705            zht(:,:) = zht(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
706         END DO
707         zwv(:,:) = 0._wp
708         DO jj = 1, jpjm1
709            DO ji = 1, fs_jpim1   ! vector opt.
710               zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e12v(ji,jj)                             &
711                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji,jj+1) * sshn(ji,jj+1) )
712            END DO
713         END DO
714         CALL lbc_lnk( zwv(:,:), 'V', 1._wp )
715         z_tmax = MAXVAL( vmask(:,:,1) * vmask_i(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) )
716         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
717         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(fse3v_n))) =', z_tmax
718         !
719         zht(:,:) = 0.0_wp
720         DO jk = 1, jpkm1
721            DO jj = 1, jpjm1
722               zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
723            END DO
724         END DO
725         zwu(:,:) = 0._wp
726         DO jj = 1, jpjm1
727            DO ji = 1, fs_jpim1   ! vector opt.
728               zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e12f(ji,jj)  &
729                        &          * (  e12t(ji  ,jj) * sshn(ji  ,jj) + e12t(ji  ,jj+1) * sshn(ji  ,jj+1) &
730                        &             + e12t(ji+1,jj) * sshn(ji+1,jj) + e12t(ji+1,jj+1) * sshn(ji+1,jj+1) )
731            END DO
732         END DO
733         CALL lbc_lnk( zht(:,:), 'F', 1._wp )
734         CALL lbc_lnk( zwu(:,:), 'F', 1._wp )
735         z_tmax = MAXVAL( fmask(:,:,1) * fmask_i(:,:) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) )
736         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
737         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(fse3f_n))) =', z_tmax
738         !
739      END IF
740
741      ! *********************************** !
742      ! After scale factors at u- v- points !
743      ! *********************************** !
744
745      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
746      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
747
748      ! *********************************** !
749      ! After depths at u- v points         !
750      ! *********************************** !
751
752      hu_a(:,:) = 0._wp                        ! Ocean depth at U-points
753      hv_a(:,:) = 0._wp                        ! Ocean depth at V-points
754      DO jk = 1, jpkm1
755         hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
756         hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
757      END DO
758      !                                        ! Inverse of the local depth
759      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
760      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
761
762      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
763      CALL wrk_dealloc( jpi, jpj, jpk, ze3t, ztu, ztv           )
764
765      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
766
767   END SUBROUTINE dom_vvl_sf_nxt
768
769
770   SUBROUTINE dom_vvl_sf_swp( kt )
771      !!----------------------------------------------------------------------
772      !!                ***  ROUTINE dom_vvl_sf_swp  ***
773      !!                   
774      !! ** Purpose :  compute time filter and swap of scale factors
775      !!               compute all depths and related variables for next time step
776      !!               write outputs and restart file
777      !!
778      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
779      !!               - reconstruct scale factor at other grid points (interpolate)
780      !!               - recompute depths and water height fields
781      !!
782      !! ** Action  :  - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
783      !!               - Recompute:
784      !!                    fse3(u/v)_b       
785      !!                    fse3w_n           
786      !!                    fse3(u/v)w_b     
787      !!                    fse3(u/v)w_n     
788      !!                    fsdept_n, fsdepw_n  and fsde3w_n
789      !!                    h(u/v) and h(u/v)r
790      !!
791      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
792      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
793      !!----------------------------------------------------------------------
794      !! * Arguments
795      INTEGER, INTENT( in )               :: kt       ! time step
796      !! * Local declarations
797      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwrk
798      INTEGER                             :: jk       ! dummy loop indices
799      !!----------------------------------------------------------------------
800
801      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
802      !
803      IF( kt == nit000 )   THEN
804         IF(lwp) WRITE(numout,*)
805         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
806         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
807      ENDIF
808      !
809      ! Time filter and swap of scale factors
810      ! =====================================
811      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
812      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
813         IF( neuler == 0 .AND. kt == nit000 ) THEN
814            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
815         ELSE
816            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
817            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
818         ENDIF
819         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
820      ENDIF
821
822      fsdept_b(:,:,:) = fsdept_n(:,:,:)
823      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
824
825      fse3t_n(:,:,:) = fse3t_a(:,:,:)
826      fse3u_n(:,:,:) = fse3u_a(:,:,:)
827      fse3v_n(:,:,:) = fse3v_a(:,:,:)
828
829      ! Compute all missing vertical scale factor and depths
830      ! ====================================================
831      ! Horizontal scale factor interpolations
832      ! --------------------------------------
833      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
834      ! - JC - hu_b, hv_b, hur_b, hvr_b also
835      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n (:,:,:), 'F'  )
836      ! Vertical scale factor interpolations
837      ! ------------------------------------
838      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
839      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
840      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
841      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
842      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
843      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
844      ! t- and w- points depth
845      ! ----------------------
846      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
847      fsdepw_n(:,:,1) = 0.0_wp
848      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
849      DO jk = 2, jpk
850         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
851         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
852         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
853      END DO
854      ! Local depth and Inverse of the local depth of the water column at u- and v- points
855      ! ----------------------------------------------------------------------------------
856      hu (:,:) = hu_a (:,:)
857      hv (:,:) = hv_a (:,:)
858
859      ! Inverse of the local depth
860      hur(:,:) = hur_a(:,:)
861      hvr(:,:) = hvr_a(:,:)
862
863      ! Local depth of the water column at t- points
864      ! --------------------------------------------
865      ht(:,:) = 0.
866      DO jk = 1, jpkm1
867         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
868      END DO
869
870      ! Write additional diagnostics
871      ! ============================
872      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) CALL dom_vvl_dia( kt) 
873
874      ! write restart file
875      ! ==================
876      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
877      !
878      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
879
880   END SUBROUTINE dom_vvl_sf_swp
881
882   SUBROUTINE dom_vvl_dia( kt )
883      !!----------------------------------------------------------------------
884      !!                ***  ROUTINE dom_vvl_dia  ***
885      !!                   
886      !! ** Purpose :  Output some diagnostics in ztilde/zlayer cases
887      !!
888      !!----------------------------------------------------------------------
889      !! * Arguments
890      INTEGER, INTENT( in )               :: kt       ! time step
891      !! * Local declarations
892      INTEGER                             :: ji,jj,jk ! dummy loop indices
893      REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, ztmp1, z2dt
894      REAL(wp), DIMENSION(4)              :: zr1
895      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zout
896      !!----------------------------------------------------------------------
897      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_dia')
898      !
899      CALL wrk_alloc( jpi, jpj, jpk, zwdw, zout )
900      !
901      ! Compute internal interfaces depths:
902      !------------------------------------
903      IF ( iom_use("dh_tilde").OR.iom_use("depw_tilde").OR.iom_use("stiff_tilde")) THEN
904         zwdw(:,:,1) = 0.e0
905         DO jj = 1, jpj
906            DO ji = 1, jpi
907               DO jk = 2, jpkm1
908                  zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + & 
909                                 & (tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
910               END DO
911            END DO
912         END DO
913      ENDIF
914      !
915      ! Output interface depth anomaly:
916      ! -------------------------------
917      IF ( iom_use("depw_tilde") ) CALL iom_put( "depw_tilde", (zwdw(:,:,:)-gdepw_0(:,:,:))*tmask(:,:,:) )
918      !
919      ! Output grid stiffness (T-points):
920      ! ---------------------------------
921      IF ( iom_use("stiff_tilde"  ) ) THEN
922         zr1(:)   = 0.e0
923         zout(:,:,jpk) = 0.e0     
924         DO ji = 2, jpim1
925            DO jj = 2, jpjm1
926               DO jk = 1, jpkm1
927                  zr1(1) = umask(ji-1,jj  ,jk) *abs( (zwdw(ji  ,jj  ,jk  )-zwdw(ji-1,jj  ,jk  )  & 
928                       &                             +zwdw(ji  ,jj  ,jk+1)-zwdw(ji-1,jj  ,jk+1)) &
929                       &                            /(zwdw(ji  ,jj  ,jk  )+zwdw(ji-1,jj  ,jk  )  &
930                       &                             -zwdw(ji  ,jj  ,jk+1)-zwdw(ji-1,jj  ,jk+1) + rsmall) )
931                  zr1(2) = umask(ji  ,jj  ,jk) *abs( (zwdw(ji+1,jj  ,jk  )-zwdw(ji  ,jj  ,jk  )  &
932                       &                             +zwdw(ji+1,jj  ,jk+1)-zwdw(ji  ,jj  ,jk+1)) &
933                       &                            /(zwdw(ji+1,jj  ,jk  )+zwdw(ji  ,jj  ,jk  )  &
934                       &                             -zwdw(ji+1,jj  ,jk+1)-zwdw(ji  ,jj  ,jk+1) + rsmall) )
935                  zr1(3) = vmask(ji  ,jj  ,jk) *abs( (zwdw(ji  ,jj+1,jk  )-zwdw(ji  ,jj  ,jk  )  &
936                       &                             +zwdw(ji  ,jj+1,jk+1)-zwdw(ji  ,jj  ,jk+1)) &
937                       &                            /(zwdw(ji  ,jj+1,jk  )+zwdw(ji  ,jj  ,jk  )  &
938                       &                             -zwdw(ji  ,jj+1,jk+1)-zwdw(ji  ,jj  ,jk+1) + rsmall) )
939                  zr1(4) = vmask(ji  ,jj-1,jk) *abs( (zwdw(ji  ,jj  ,jk  )-zwdw(ji  ,jj-1,jk  )  &
940                       &                             +zwdw(ji  ,jj  ,jk+1)-zwdw(ji  ,jj-1,jk+1)) &
941                       &                            /(zwdw(ji  ,jj  ,jk  )+zwdw(ji  ,jj-1,jk  )  &
942                       &                             -zwdw(ji,  jj  ,jk+1)-zwdw(ji  ,jj-1,jk+1) + rsmall) )
943                  zout(ji,jj,jk) = MAXVAL(zr1(1:4))
944               END DO
945            END DO
946         END DO
947
948         CALL lbc_lnk( zout, 'T', 1. )
949         CALL iom_put( "stiff_tilde", zout(:,:,:) ) 
950      END IF
951      ! Output Horizontal Laplacian of interfaces depths (W-points):
952      ! ------------------------------------------------------------
953      IF ( iom_use("dh_tilde")   ) THEN
954         !
955         zout(:,:,1  )=0._wp
956         zout(:,:,jpk)=0._wp
957         DO jk = 2, jpkm1
958            DO jj = 2, jpjm1
959               DO ji = fs_2, fs_jpim1   ! vector opt.
960                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) )
961                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) )
962                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) )
963                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) )
964!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj)) 
965                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 
966!                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)
967                  zout(ji,jj,jk) = ztmp1
968               END DO
969            END DO           
970         END DO
971         ! Mask open boundaries:
972#if defined key_bdy
973         IF (lk_bdy) THEN
974            DO jk = 1, jpkm1
975               zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:)
976            END DO
977         ENDIF
978#endif
979         CALL lbc_lnk( zout(:,:,:), 'T', 1. )
980         zwdw(:,:,:) = zout(:,:,:)
981         DO jk = 2, jpkm1
982            DO jj = 2, jpjm1
983               DO ji = fs_2, fs_jpim1   ! vector opt.
984                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) )
985                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) )
986                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) )
987                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) )
988!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj)) 
989                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 
990                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)
991               END DO
992            END DO           
993         END DO
994         ! Mask open boundaries:
995#if defined key_bdy
996         IF (lk_bdy) THEN
997            DO jk = 1, jpkm1
998               zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:)
999            END DO
1000         ENDIF
1001#endif
1002         CALL lbc_lnk( zout(:,:,:), 'T', 1. )
1003         !
1004         CALL iom_put( "dh_tilde", zout(:,:,:) )
1005      ENDIF
1006      !
1007      ! Output vertical Laplacian of interfaces depths (W-points):
1008      ! ----------------------------------------------------------
1009      IF ( iom_use("dz_tilde"  ) ) THEN
1010         zout(:,:,1  ) = 0._wp
1011         zout(:,:,jpk) = 0._wp
1012         DO jk=2,jpkm1
1013            zout(:,:,jk) = 2._wp*ABS(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)-tilde_e3t_n(:,:,jk-1)-e3t_0(:,:,jk-1)) & 
1014                               &   /(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)+tilde_e3t_n(:,:,jk-1)+e3t_0(:,:,jk-1)) &
1015                               &   * tmask(:,:,jk)
1016         END DO
1017         CALL iom_put( "dz_tilde", zout(:,:,:) ) 
1018      END IF
1019      !
1020      !
1021      ! Output low pass U-velocity:
1022      ! ---------------------------
1023      IF ( iom_use("un_lf_tilde"  ).AND.ln_vvl_ztilde ) THEN
1024         zout(:,:,jpk) = 0.e0 
1025         DO jk=1,jpkm1
1026            zout(:,:,jk) = un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:)
1027         END DO
1028         CALL iom_put( "un_lf_tilde", zout(:,:,:) )
1029      END IF
1030      !
1031      ! Output low pass V-velocity:
1032      ! ---------------------------
1033      IF ( iom_use("vn_lf_tilde"  ).AND.ln_vvl_ztilde ) THEN
1034         zout(:,:,jpk) = 0.e0 
1035         DO jk=1,jpkm1
1036            zout(:,:,jk) = vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:)
1037         END DO
1038         CALL iom_put( "vn_lf_tilde", zout(:,:,:) )
1039      END IF
1040      !
1041      CALL wrk_dealloc( jpi, jpj, jpk, zwdw, zout )
1042      !
1043      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_dia')
1044      !
1045   END SUBROUTINE dom_vvl_dia
1046
1047   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
1048      !!---------------------------------------------------------------------
1049      !!                  ***  ROUTINE dom_vvl__interpol  ***
1050      !!
1051      !! ** Purpose :   interpolate scale factors from one grid point to another
1052      !!
1053      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
1054      !!                - horizontal interpolation: grid cell surface averaging
1055      !!                - vertical interpolation: simple averaging
1056      !!----------------------------------------------------------------------
1057      !! * Arguments
1058      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
1059      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
1060      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
1061      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1062      !! * Local declarations
1063      INTEGER ::   nmet                                                ! horizontal interpolation method
1064      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
1065      INTEGER ::   jkbot                                               ! bottom level
1066      LOGICAL ::   l_is_orca                                           ! local logical
1067      REAL(wp) :: zmin, zdo, zup, ztap
1068      REAL(wp), POINTER, DIMENSION(:,:)   :: zs                        ! Surface interface depth anomaly
1069      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw                        ! Interface depth anomaly
1070      !!----------------------------------------------------------------------
1071      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
1072         !
1073      l_is_orca = .FALSE.
1074      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
1075
1076 
1077      nmet=1 ! Original method (Surely wrong)
1078!      nmet=1 ! Interface interpolation
1079!      nmet=2 ! Internal interfaces interpolation only, spread barotropic increment
1080!      nmet=1
1081      ztap = 0.1_wp ! Minimum fraction of T-point thickness at cell interfaces
1082
1083      IF ( (nmet==1).OR.(nmet==2) ) THEN
1084         SELECT CASE ( pout )
1085            !
1086         CASE( 'U', 'V', 'F' )
1087            ! Compute interface depth anomaly at T-points
1088            CALL wrk_alloc( jpi, jpj, jpk, zw )
1089            CALL wrk_alloc( jpi, jpj, zs )
1090            !
1091            zw(:,:,:) =  0._wp   
1092            !
1093            DO jj = 1, jpj
1094               DO ji = 1, jpi
1095                  jkbot = mbkt(ji,jj)
1096                  DO jk=jkbot,1,-1
1097                     zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk)
1098                  END DO
1099               END DO
1100            END DO 
1101            !
1102            !
1103            IF (nmet==2) THEN        ! Consider "internal" interfaces only
1104               zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh)
1105               !
1106               zw(:,:,:) = 0._wp
1107               DO jj = 1, jpj
1108                  DO ji = 1, jpi
1109                     jkbot = mbkt(ji,jj)
1110                     DO jk=jkbot,1,-1
1111                        zw(ji,jj,jk) = zw(ji,jj,jk+1) - ( pe3_in(ji,jj,jk)                     &
1112                                     & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1))  &
1113                                     & - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
1114                     END DO
1115                  END DO
1116               END DO
1117            ENDIF 
1118            !
1119         END SELECT
1120      END IF
1121
1122      pe3_out(:,:,:) = 0._wp
1123
1124      SELECT CASE ( pout )
1125         !               ! ------------------------------------- !
1126      CASE( 'U' )        ! interpolation from T-point to U-point !
1127         !               ! ------------------------------------- !     
1128         ! horizontal surface weighted interpolation
1129         IF (nmet==0) THEN
1130            DO jk = 1, jpk
1131               DO jj = 1, jpjm1
1132                  DO ji = 1, fs_jpim1   ! vector opt.
1133                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
1134                        &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
1135                        &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
1136                  END DO
1137               END DO
1138            END DO
1139         ENDIF
1140         !
1141         IF ( (nmet==1).OR.(nmet==2) ) THEN
1142            DO jj = 1, jpjm1
1143               DO ji = 1, fs_jpim1   ! vector opt.
1144                  ! Correction at last level:
1145                  jkbot = mbku(ji,jj)
1146                  pe3_out(ji,jj,jkbot) = -0.5_wp * umask(ji,jj,jkbot) * r1_e12u(ji,jj) &
1147                        &                    * (   e12t(ji  ,jj) * zw(ji  ,jj,jkbot)   &
1148                        &                        + e12t(ji+1,jj) * zw(ji+1,jj,jkbot)   )
1149                  !
1150                  ! If there is a step, taper bottom interface:
1151                  IF (hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ) THEN
1152                     IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
1153!                        zmin = ztap * pe3_in(ji+1,jj,jkbot)
1154                        zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
1155                     ELSE
1156!                        zmin = ztap * pe3_in(ji  ,jj,jkbot)
1157                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
1158                     ENDIF
1159                     zmin = -e3u_0(ji,jj,jkbot) + zmin
1160                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
1161                  ENDIF
1162                  !
1163                  zdo =  -pe3_out(ji,jj,jkbot)
1164                  DO jk=jkbot-1,1,-1
1165                     zup = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji  ,jj) &
1166                              &                     *(   e12t(ji  ,jj) * zw(ji  ,jj,jk)  &
1167                              &                         +e12t(ji+1,jj) * zw(ji+1,jj,jk)  )
1168                     pe3_out(ji,jj,jk) = zdo - zup 
1169                     zdo =  zdo - pe3_out(ji,jj,jk)
1170                  END DO
1171               END DO
1172            END DO
1173         END IF
1174         !
1175         IF (nmet==2) THEN           ! Spread sea level anomaly
1176            DO jj = 1, jpjm1
1177               DO ji = 1, fs_jpim1   ! vector opt.
1178                  DO jk=1,jpk
1179                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                 &
1180                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )             & 
1181                           &               / ( hu_0(ji,jj) + 1._wp - umask_i(ji,jj) )            &
1182                           &               * 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)           &
1183                           &               * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji+1,jj)*zs(ji+1,jj) )       
1184                  END DO
1185               END DO
1186            END DO
1187            !
1188         ENDIF
1189         !
1190         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1191         ! boundary conditions
1192         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
1193         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
1194         !
1195         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
1196         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
1197         !
1198         !               ! ------------------------------------- !
1199      CASE( 'V' )        ! interpolation from T-point to V-point !
1200         !               ! ------------------------------------- !
1201         ! horizontal surface weighted interpolation
1202         IF (nmet==0) THEN
1203            DO jk = 1, jpk
1204               DO jj = 1, jpjm1
1205                  DO ji = 1, fs_jpim1   ! vector opt.
1206                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
1207                        &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
1208                        &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
1209                  END DO
1210               END DO
1211            END DO
1212         ENDIF
1213         !
1214         IF ( (nmet==1).OR.(nmet==2) ) THEN
1215            DO jj = 1, jpjm1
1216               DO ji = 1, fs_jpim1   ! vector opt.
1217                  ! Correction at last level:
1218                  jkbot = mbkv(ji,jj)
1219                  pe3_out(ji,jj,jkbot) = -0.5_wp * vmask(ji,jj,jkbot) * r1_e12v(ji,jj) &
1220                        &                    * (   e12t(ji,jj  ) * zw(ji,jj  ,jkbot)   &
1221                        &                        + e12t(ji,jj+1) * zw(ji,jj+1,jkbot)   )
1222                  !
1223                  ! If there is a step, taper bottom interface:
1224                  IF (hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ) THEN
1225                     IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
1226!                        zmin = ztap * pe3_in(ji,jj+1,jkbot)
1227                        zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
1228                     ELSE
1229!                        zmin = ztap * pe3_in(ji,jj  ,jkbot)
1230                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
1231                     ENDIF
1232                     zmin = -e3v_0(ji,jj,jkbot) + zmin
1233                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
1234                  ENDIF
1235                  !
1236                  zdo =  -pe3_out(ji,jj,jkbot)
1237                  DO jk=jkbot-1,1,-1
1238                     zup = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj  ) &
1239                              &                     *  ( e12t(ji,jj  ) * zw(ji,jj  ,jk) &
1240                              &                         +e12t(ji,jj+1) * zw(ji,jj+1,jk) )
1241                     !
1242                     pe3_out(ji,jj,jk) = zdo - zup 
1243                     zdo =  zdo - pe3_out(ji,jj,jk)
1244                  END DO
1245               END DO
1246            END DO
1247         END IF
1248         !
1249         IF (nmet==2) THEN           ! Spread sea level anomaly
1250            !
1251            DO jj = 1, jpjm1
1252               DO ji = 1, fs_jpim1   ! vector opt.
1253                  DO jk=1,jpk
1254                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                 &
1255                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )             & 
1256                           &               / ( hv_0(ji,jj) + 1._wp - vmask_i(ji,jj) )            &
1257                           &               * 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)           &
1258                           &               * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji,jj+1)*zs(ji,jj+1) )       
1259                  END DO
1260               END DO
1261            END DO
1262            !
1263         ENDIF
1264         !
1265         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1266         ! boundary conditions
1267         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
1268         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
1269         !
1270         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
1271         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
1272         !
1273         !               ! ------------------------------------- !
1274      CASE( 'F' )        ! interpolation from T-point to F-point !
1275         !               ! ------------------------------------- !
1276         ! horizontal surface weighted interpolation
1277         IF (nmet==0) THEN
1278            DO jk=1,jpk
1279               DO jj = 1, jpjm1
1280                  DO ji = 1, fs_jpim1   ! vector opt.
1281                     pe3_out(ji,jj,jk) = 0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
1282                           &      * (  e12t(ji  ,jj  ) * ( pe3_in(ji  ,jj  ,jk)-e3t_0(ji  ,jj  ,jk) )   & 
1283                           &         + e12t(ji  ,jj+1) * ( pe3_in(ji  ,jj+1,jk)-e3t_0(ji  ,jj+1,jk) )   &
1284                           &         + e12t(ji+1,jj  ) * ( pe3_in(ji+1,jj  ,jk)-e3t_0(ji+1,jj  ,jk) )   & 
1285                           &         + e12t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ))                                                     
1286                  END DO
1287               END DO
1288            END DO
1289         ENDIF
1290         !
1291         IF ( (nmet==1).OR.(nmet==2) ) THEN
1292            DO jj = 1, jpjm1
1293               DO ji = 1, fs_jpim1   ! vector opt.
1294                  ! bottom correction:
1295                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 
1296                  pe3_out(ji,jj,jkbot) = -0.25_wp * umask(ji,jj,jkbot) * umask(ji,jj+1,jkbot) * r1_e12f(ji,jj) &
1297                        &                         * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jkbot)   &
1298                        &                            + e12t(ji+1,jj  ) * zw(ji+1,jj  ,jkbot)   &
1299                        &                            + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jkbot)   &
1300                        &                            + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jkbot)   )
1301                  !
1302                  ! If there is a step, taper bottom interface:
1303                  IF (hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ) THEN
1304                     IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
1305                        IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN
1306!                           zmin = ztap * pe3_in(ji+1,jj+1,jkbot)
1307                           zmin = ztap * (-zw(ji+1,jj+1,jkbot)+e3t_0(ji+1,jj+1,jkbot))
1308                        ELSE
1309!                           zmin = ztap * pe3_in(ji  ,jj+1,jkbot)
1310                           zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
1311                        ENDIF
1312                     ELSE
1313                        IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
1314!                           zmin = ztap * pe3_in(ji+1,jj  ,jkbot)
1315                           zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
1316                        ELSE
1317!                           zmin = ztap * pe3_in(ji  ,jj  ,jkbot)
1318                           zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
1319                        ENDIF
1320                     ENDIF
1321                     zmin = -e3f_0(ji,jj,jkbot) + zmin
1322                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
1323                  ENDIF
1324                  !
1325                  zdo =  -pe3_out(ji,jj,jkbot)
1326                  DO jk=jkbot-1,1,-1
1327                     zup =  0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
1328                           &        * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  & 
1329                           &           + e12t(ji+1,jj  ) * zw(ji+1,jj  ,jk)  & 
1330                           &           + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  & 
1331                           &           + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  )
1332                     pe3_out(ji,jj,jk) = zdo - zup 
1333                     !
1334                     zdo =  zdo - pe3_out(ji,jj,jk)                                   
1335                  END DO
1336                  !
1337               END DO
1338            END DO
1339         ENDIF
1340         !
1341         IF (nmet==2) THEN           ! Spread sea level anomaly
1342            !
1343            DO jj = 1, jpjm1
1344               DO ji = 1, fs_jpim1   ! vector opt.
1345                  DO jk=1,jpk
1346                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                         &
1347                           &               + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) )                     & 
1348                           &               / ( hf_0(ji,jj) + 1._wp - umask_i(ji,jj)*umask_i(ji,jj+1) )   &
1349                           &               * 0.25_wp * umask(ji,jj,jk)*umask(ji,jj+1,jk)*r1_e12f(ji,jj)  &
1350                           &               * ( e12t(ji  ,jj)*zs(ji  ,jj) + e12t(ji  ,jj+1)*zs(ji  ,jj+1) &
1351                           &                  +e12t(ji+1,jj)*zs(ji+1,jj) + e12t(ji+1,jj+1)*zs(ji+1,jj+1) )               
1352                  END DO
1353               END DO
1354            END DO
1355         END IF
1356         !
1357         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1358         ! boundary conditions
1359         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
1360         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
1361         !
1362         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
1363         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
1364         !
1365         !               ! ------------------------------------- !
1366      CASE( 'W' )        ! interpolation from T-point to W-point !
1367         !               ! ------------------------------------- !
1368         ! vertical simple interpolation
1369         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
1370         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1371         DO jk = 2, jpk
1372            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
1373               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
1374         END DO
1375         !               ! -------------------------------------- !
1376      CASE( 'UW' )       ! interpolation from U-point to UW-point !
1377         !               ! -------------------------------------- !
1378         ! vertical simple interpolation
1379         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
1380         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1381         DO jk = 2, jpk
1382            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
1383               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
1384         END DO
1385         !               ! -------------------------------------- !
1386      CASE( 'VW' )       ! interpolation from V-point to VW-point !
1387         !               ! -------------------------------------- !
1388         ! vertical simple interpolation
1389         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
1390         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1391         DO jk = 2, jpk
1392            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
1393               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
1394         END DO
1395      END SELECT
1396      !
1397
1398      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
1399
1400   END SUBROUTINE dom_vvl_interpol
1401
1402   SUBROUTINE dom_vvl_rst( kt, cdrw )
1403      !!---------------------------------------------------------------------
1404      !!                   ***  ROUTINE dom_vvl_rst  ***
1405      !!                     
1406      !! ** Purpose :   Read or write VVL file in restart file
1407      !!
1408      !! ** Method  :   use of IOM library
1409      !!                if the restart does not contain vertical scale factors,
1410      !!                they are set to the _0 values
1411      !!                if the restart does not contain vertical scale factors increments (z_tilde),
1412      !!                they are set to 0.
1413      !!----------------------------------------------------------------------
1414      !! * Arguments
1415      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1416      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1417      !! * Local declarations
1418      INTEGER ::   jk
1419      INTEGER, DIMENSION(7) ::   id            ! local integers
1420      REAL(wp), POINTER, DIMENSION(:,:) :: zhdiv
1421      !!----------------------------------------------------------------------
1422      !
1423      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
1424      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1425         !                                   ! ===============
1426         IF( ln_rstart ) THEN                   !* Read the restart file
1427            CALL rst_read_open                  !  open the restart file if necessary
1428            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
1429            !
1430            id(1) = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
1431            id(2) = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
1432            id(3) = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
1433            id(4) = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
1434            id(5) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
1435            id(6) = iom_varid( numror, 'un_lf'   , ldstop = .FALSE. )
1436            id(7) = iom_varid( numror, 'vn_lf'   , ldstop = .FALSE. )
1437
1438            !                             ! --------- !
1439            !                             ! all cases !
1440            !                             ! --------- !
1441            IF( MIN( id(1), id(2) ) > 0 ) THEN       ! all required arrays exist
1442               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
1443               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
1444               ! needed to restart if land processor not computed
1445               IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
1446               WHERE ( tmask(:,:,:) == 0.0_wp ) 
1447                  fse3t_n(:,:,:) = e3t_0(:,:,:)
1448                  fse3t_b(:,:,:) = e3t_0(:,:,:)
1449               END WHERE
1450               IF( neuler == 0 ) THEN
1451                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
1452               ENDIF
1453            ELSE IF( id(1) > 0 ) THEN
1454               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
1455               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
1456               IF(lwp) write(numout,*) 'neuler is forced to 0'
1457               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
1458               fse3t_b(:,:,:) = fse3t_n(:,:,:)
1459               neuler = 0
1460            ELSE IF( id(2) > 0 ) THEN
1461               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
1462               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
1463               IF(lwp) write(numout,*) 'neuler is forced to 0'
1464               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
1465               fse3t_b(:,:,:) = fse3t_n(:,:,:)
1466               neuler = 0
1467            ELSE
1468               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
1469               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1470               IF(lwp) write(numout,*) 'neuler is forced to 0'
1471               DO jk=1,jpk
1472                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
1473                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
1474                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
1475               END DO
1476               fse3t_b(:,:,:) = fse3t_n(:,:,:)
1477               neuler = 0
1478            ENDIF
1479            !                             ! ----------- !
1480            IF( ln_vvl_zstar ) THEN       ! z_star case !
1481               !                          ! ----------- !
1482               IF( MIN( id(3), id(4) ) > 0 ) THEN
1483                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
1484               ENDIF
1485               !                          ! ----------------------- !
1486            ELSE                          ! z_tilde and layer cases !
1487               !                          ! ----------------------- !
1488               IF( MIN( id(3), id(4) ) > 0 ) THEN  ! all required arrays exist
1489                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
1490                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
1491               ELSE                            ! one at least array is missing
1492                  tilde_e3t_b(:,:,:) = 0.0_wp
1493                  tilde_e3t_n(:,:,:) = 0.0_wp
1494                  !
1495                  neuler = 0
1496               ENDIF                      ! ------------ !
1497               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
1498                  !                       ! ------------ !
1499                  IF( MINVAL(id(5:7) ) > 0 ) THEN  ! all required arrays exist
1500                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) )
1501                     CALL iom_get( numror, jpdom_autoglo, 'un_lf'   ,    un_lf(:,:,:) )
1502                     CALL iom_get( numror, jpdom_autoglo, 'vn_lf'   ,    vn_lf(:,:,:) )
1503                  ELSE                            ! one at least array is missing
1504                     hdivn_lf(:,:,:) = 0.0_wp
1505                     un_lf(:,:,:)    = 0.0_wp
1506                     vn_lf(:,:,:)    = 0.0_wp
1507                     neuler = 0
1508                  ENDIF
1509               ENDIF
1510               !
1511            ENDIF
1512            !
1513         ELSE                                   !* Initialize at "rest"
1514            fse3t_b(:,:,:) = e3t_0(:,:,:)
1515            fse3t_n(:,:,:) = e3t_0(:,:,:)
1516            sshn(:,:) = 0.0_wp
1517            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1518               tilde_e3t_b(:,:,:) = 0.0_wp
1519               tilde_e3t_n(:,:,:) = 0.0_wp
1520            END IF
1521            IF( ln_vvl_ztilde ) THEN
1522               hdivn_lf(:,:,:) = 0.0_wp
1523               un_lf(:,:,:)    = 0.0_wp
1524               vn_lf(:,:,:)    = 0.0_wp
1525            ENDIF
1526         ENDIF
1527
1528      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1529         !                                   ! ===================
1530         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1531         !                                           ! --------- !
1532         !                                           ! all cases !
1533         !                                           ! --------- !
1534         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
1535         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
1536         !                                           ! ----------------------- !
1537         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
1538            !                                        ! ----------------------- !
1539            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
1540            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
1541         END IF
1542         !                                           ! -------------!   
1543         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
1544            !                                        ! ------------ !
1545            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:) )
1546            CALL iom_rstput( kt, nitrst, numrow, 'un_lf'   , un_lf(:,:,:)    )
1547            CALL iom_rstput( kt, nitrst, numrow, 'vn_lf'   , vn_lf(:,:,:)    )
1548         ENDIF
1549
1550      ENDIF
1551      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
1552
1553   END SUBROUTINE dom_vvl_rst
1554
1555   SUBROUTINE dom_vvl_ctl
1556      !!---------------------------------------------------------------------
1557      !!                  ***  ROUTINE dom_vvl_ctl  ***
1558      !!               
1559      !! ** Purpose :   Control the consistency between namelist options
1560      !!                for vertical coordinate
1561      !!----------------------------------------------------------------------
1562      INTEGER ::   ioptio
1563      INTEGER ::   ios
1564
1565      NAMELIST/nam_vvl/ ln_vvl_zstar               , ln_vvl_ztilde                       , &
1566                      & ln_vvl_layer               , ln_vvl_ztilde_as_zstar              , &
1567                      & ln_vvl_zstar_at_eqtor      , ln_vvl_zstar_on_shelf               , &
1568                      & ln_vvl_adv_cn2             , ln_vvl_adv_fct                      , &
1569                      & ln_vvl_lap                 , ln_vvl_blp                          , &
1570                      & rn_ahe3_lap                , rn_ahe3_blp                         , &
1571                      & rn_rst_e3t                 , rn_lf_cutoff                        , &
1572                      & ln_vvl_regrid              , rn_zdef_max                         , &
1573                      & ln_vvl_dbg   ! not yet implemented: ln_vvl_kepe
1574      !!----------------------------------------------------------------------
1575
1576      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
1577      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1578901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
1579
1580      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
1581      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1582902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
1583      IF(lwm) WRITE ( numond, nam_vvl )
1584
1585      IF(lwp) THEN                    ! Namelist print
1586         WRITE(numout,*)
1587         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1588         WRITE(numout,*) '~~~~~~~~~~~'
1589         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
1590         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1591         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1592         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
1593         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1594         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
1595         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
1596         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1597         WRITE(numout,*) '      ztilde on shelves          ln_vvl_zstar_on_shelf  = ', ln_vvl_zstar_on_shelf
1598         WRITE(numout,*) '           Namelist nam_vvl : thickness advection scheme'
1599         WRITE(numout,*) '              2nd order                  ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2
1600         WRITE(numout,*) '              2nd order FCT              ln_vvl_adv_fct = ', ln_vvl_adv_fct                   
1601         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion scheme'
1602         WRITE(numout,*) '              Laplacian                  ln_vvl_lap     = ', ln_vvl_lap
1603         WRITE(numout,*) '              Bilaplacian                ln_vvl_blp     = ', ln_vvl_blp
1604         WRITE(numout,*) '              Laplacian   coefficient    rn_ahe3_lap    = ', rn_ahe3_lap
1605         WRITE(numout,*) '              Bilaplacian coefficient    rn_ahe3_blp    = ', rn_ahe3_blp
1606         WRITE(numout,*) '           Namelist nam_vvl : layers regriding'
1607         WRITE(numout,*) '              ln_vvl_regrid                              = ', ln_vvl_regrid
1608         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
1609         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
1610         IF( ln_vvl_ztilde_as_zstar ) THEN
1611            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
1612            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
1613            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
1614            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
1615            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1616            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
1617         ELSE
1618            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
1619            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
1620            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
1621            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
1622         ENDIF
1623         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
1624         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
1625      ENDIF
1626
1627      ioptio = 0                      ! Parameter control
1628      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1629      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
1630      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
1631      IF( ln_vvl_layer           )        ioptio = ioptio + 1
1632
1633      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1634
1635      IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN
1636         ioptio = 0                      ! Choose one advection scheme at most
1637         IF( ln_vvl_adv_cn2         )        ioptio = ioptio + 1
1638         IF( ln_vvl_adv_fct         )        ioptio = ioptio + 1
1639         IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' )
1640      ENDIF
1641
1642      IF(lwp) THEN                   ! Print the choice
1643         WRITE(numout,*)
1644         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
1645         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
1646         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
1647         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
1648         ! - ML - Option not developed yet
1649         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
1650         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
1651      ENDIF
1652
1653#if defined key_agrif
1654      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
1655#endif
1656
1657   END SUBROUTINE dom_vvl_ctl
1658
1659   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1660      !!---------------------------------------------------------------------
1661      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
1662      !!                     
1663      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
1664      !!                scale factors at locations that have been individually
1665      !!                modified in domhgr. Such modifications break the
1666      !!                relationship between e12t and e1u*e2u etc.
1667      !!                Recompute some scale factors ignoring the modified metric.
1668      !!----------------------------------------------------------------------
1669      !! * Arguments
1670      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
1671      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
1672      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
1673      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1674      !! * Local declarations
1675      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
1676      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
1677      INTEGER ::   isrow                                               ! index for ORCA1 starting row
1678      !! acc
1679      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1680      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1681      !!
1682      !                                                ! =====================
1683      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
1684         !                                             ! =====================
1685      !! acc
1686         IF( nn_cla == 0 ) THEN
1687            !
1688            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
1689            ij0 = 102   ;   ij1 = 102
1690            DO jk = 1, jpkm1
1691               DO jj = mj0(ij0), mj1(ij1)
1692                  DO ji = mi0(ii0), mi1(ii1)
1693                     SELECT CASE ( pout )
1694                     CASE( 'U' )
1695                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1696                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1697                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1698                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1699                     CASE( 'F' )
1700                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1701                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1702                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1703                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1704                     END SELECT
1705                  END DO
1706               END DO
1707            END DO
1708            !
1709            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
1710            ij0 =  88   ;   ij1 =  88
1711            DO jk = 1, jpkm1
1712               DO jj = mj0(ij0), mj1(ij1)
1713                  DO ji = mi0(ii0), mi1(ii1)
1714                     SELECT CASE ( pout )
1715                     CASE( 'U' )
1716                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1717                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1718                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1719                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1720                     CASE( 'V' )
1721                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1722                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1723                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1724                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1725                     CASE( 'F' )
1726                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1727                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1728                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1729                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1730                     END SELECT
1731                  END DO
1732               END DO
1733            END DO
1734         ENDIF
1735
1736         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
1737         ij0 = 116   ;   ij1 = 116
1738         DO jk = 1, jpkm1
1739            DO jj = mj0(ij0), mj1(ij1)
1740               DO ji = mi0(ii0), mi1(ii1)
1741                  SELECT CASE ( pout )
1742                  CASE( 'U' )
1743                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1744                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1745                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1746                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1747                  CASE( 'F' )
1748                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1749                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1750                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1751                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1752                  END SELECT
1753               END DO
1754            END DO
1755         END DO
1756      ENDIF
1757      !
1758         !                                             ! =====================
1759      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1760         !                                             ! =====================
1761         ! This dirty section will be suppressed by simplification process:
1762         ! all this will come back in input files
1763         ! Currently these hard-wired indices relate to configuration with
1764         ! extend grid (jpjglo=332)
1765         ! which had a grid-size of 362x292.
1766         isrow = 332 - jpjglo
1767         !
1768         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified)
1769         ij0 = 241 - isrow   ;   ij1 = 241 - isrow
1770         DO jk = 1, jpkm1
1771            DO jj = mj0(ij0), mj1(ij1)
1772               DO ji = mi0(ii0), mi1(ii1)
1773                  SELECT CASE ( pout )
1774                  CASE( 'U' )
1775                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1776                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1777                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1778                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1779                  CASE( 'F' )
1780                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1781                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1782                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1783                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1784                  END SELECT
1785               END DO
1786            END DO
1787         END DO
1788         !
1789         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1790         ij0 = 248 - isrow   ;   ij1 = 248 - isrow
1791         DO jk = 1, jpkm1
1792            DO jj = mj0(ij0), mj1(ij1)
1793               DO ji = mi0(ii0), mi1(ii1)
1794                  SELECT CASE ( pout )
1795                  CASE( 'U' )
1796                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1797                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1798                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1799                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1800                  CASE( 'F' )
1801                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1802                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1803                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1804                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1805                  END SELECT
1806               END DO
1807            END DO
1808         END DO
1809         !
1810         ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1811         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
1812         DO jk = 1, jpkm1
1813            DO jj = mj0(ij0), mj1(ij1)
1814               DO ji = mi0(ii0), mi1(ii1)
1815                  SELECT CASE ( pout )
1816                  CASE( 'V' )
1817                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1818                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1819                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1820                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1821                  END SELECT
1822               END DO
1823            END DO
1824         END DO
1825         !
1826         ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1827         ij0 = 164 - isrow   ;   ij1 = 165 - isrow
1828         DO jk = 1, jpkm1
1829            DO jj = mj0(ij0), mj1(ij1)
1830               DO ji = mi0(ii0), mi1(ii1)
1831                  SELECT CASE ( pout )
1832                  CASE( 'V' )
1833                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1834                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1835                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1836                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1837                  END SELECT
1838               END DO
1839            END DO
1840         END DO
1841         !
1842         ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1843         ij0 = 164 - isrow  ;   ij1 = 165  - isrow 
1844         DO jk = 1, jpkm1
1845            DO jj = mj0(ij0), mj1(ij1)
1846               DO ji = mi0(ii0), mi1(ii1)
1847                  SELECT CASE ( pout )
1848                  CASE( 'V' )
1849                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1850                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1851                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1852                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1853                  END SELECT
1854               END DO
1855            END DO
1856         END DO
1857         !
1858         ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified)
1859         ij0 = 164 - isrow    ;   ij1 = 165  - isrow 
1860         DO jk = 1, jpkm1
1861            DO jj = mj0(ij0), mj1(ij1)
1862               DO ji = mi0(ii0), mi1(ii1)
1863                  SELECT CASE ( pout )
1864                  CASE( 'V' )
1865                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1866                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1867                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1868                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1869                  END SELECT
1870               END DO
1871            END DO
1872         END DO
1873         !
1874         ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1875         ij0 = 181 - isrow    ;   ij1 = 182 - isrow 
1876         DO jk = 1, jpkm1
1877            DO jj = mj0(ij0), mj1(ij1)
1878               DO ji = mi0(ii0), mi1(ii1)
1879                  SELECT CASE ( pout )
1880                  CASE( 'V' )
1881                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1882                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1883                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1884                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1885                  END SELECT
1886               END DO
1887            END DO
1888         END DO
1889         !
1890         ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1891         ij0 = 181 - isrow    ;   ij1 = 182 - isrow 
1892         DO jk = 1, jpkm1
1893            DO jj = mj0(ij0), mj1(ij1)
1894               DO ji = mi0(ii0), mi1(ii1)
1895                  SELECT CASE ( pout )
1896                  CASE( 'V' )
1897                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1898                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1899                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1900                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1901                  END SELECT
1902               END DO
1903            END DO
1904         END DO
1905      ENDIF
1906         !                                             ! =====================
1907      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1908         !                                             ! =====================
1909         !
1910         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1911         ij0 = 327   ;   ij1 = 327
1912         DO jk = 1, jpkm1
1913            DO jj = mj0(ij0), mj1(ij1)
1914               DO ji = mi0(ii0), mi1(ii1)
1915                  SELECT CASE ( pout )
1916                  CASE( 'U' )
1917                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1918                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1919                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1920                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1921                  CASE( 'F' )
1922                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1923                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1924                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1925                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1926                  END SELECT
1927               END DO
1928            END DO
1929         END DO
1930         !
1931         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1932         ij0 = 343   ;   ij1 = 343
1933         DO jk = 1, jpkm1
1934            DO jj = mj0(ij0), mj1(ij1)
1935               DO ji = mi0(ii0), mi1(ii1)
1936                  SELECT CASE ( pout )
1937                  CASE( 'U' )
1938                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1939                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1940                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1941                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1942                  CASE( 'F' )
1943                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1944                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1945                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1946                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1947                  END SELECT
1948               END DO
1949            END DO
1950         END DO
1951         !
1952         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1953         ij0 = 232   ;   ij1 = 232
1954         DO jk = 1, jpkm1
1955            DO jj = mj0(ij0), mj1(ij1)
1956               DO ji = mi0(ii0), mi1(ii1)
1957                  SELECT CASE ( pout )
1958                  CASE( 'U' )
1959                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1960                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1961                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1962                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1963                  CASE( 'F' )
1964                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1965                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1966                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1967                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1968                  END SELECT
1969               END DO
1970            END DO
1971         END DO
1972         !
1973         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1974         ij0 = 232   ;   ij1 = 232
1975         DO jk = 1, jpkm1
1976            DO jj = mj0(ij0), mj1(ij1)
1977               DO ji = mi0(ii0), mi1(ii1)
1978                  SELECT CASE ( pout )
1979                  CASE( 'U' )
1980                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1981                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1982                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1983                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1984                  CASE( 'F' )
1985                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1986                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1987                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1988                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1989                  END SELECT
1990               END DO
1991            END DO
1992         END DO
1993         !
1994         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1995         ij0 = 270   ;   ij1 = 270
1996         DO jk = 1, jpkm1
1997            DO jj = mj0(ij0), mj1(ij1)
1998               DO ji = mi0(ii0), mi1(ii1)
1999                  SELECT CASE ( pout )
2000                  CASE( 'U' )
2001                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
2002                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
2003                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
2004                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
2005                  CASE( 'F' )
2006                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
2007                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
2008                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
2009                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
2010                  END SELECT
2011               END DO
2012            END DO
2013         END DO
2014         !
2015         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
2016         ij0 = 232   ;   ij1 = 233
2017         DO jk = 1, jpkm1
2018            DO jj = mj0(ij0), mj1(ij1)
2019               DO ji = mi0(ii0), mi1(ii1)
2020                  SELECT CASE ( pout )
2021                  CASE( 'V' )
2022                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
2023                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
2024                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
2025                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
2026                  END SELECT
2027               END DO
2028            END DO
2029         END DO
2030         !
2031         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
2032         ij0 = 276   ;   ij1 = 276
2033         DO jk = 1, jpkm1
2034            DO jj = mj0(ij0), mj1(ij1)
2035               DO ji = mi0(ii0), mi1(ii1)
2036                  SELECT CASE ( pout )
2037                  CASE( 'V' )
2038                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
2039                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
2040                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
2041                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
2042                  END SELECT
2043               END DO
2044            END DO
2045         END DO
2046      ENDIF
2047   END SUBROUTINE dom_vvl_orca_fix
2048
2049   SUBROUTINE dom_vvl_zdf( kt, p2dt ) 
2050      !!----------------------------------------------------------------------
2051      !!                  ***  ROUTINE dom_vvl_zdf  ***
2052      !!
2053      !! ** Purpose :   Do vertical thicknesses anomaly diffusion
2054      !!
2055      !! ** Method  : 
2056      !!
2057      !! ** Action  :
2058      !!---------------------------------------------------------------------
2059      USE oce     , ONLY:   zwd => ua       , zws => va         ! (ua,va) used as 3D workspace
2060      !
2061      INTEGER                         , INTENT(in) ::   kt     ! ocean time-step index
2062      REAL(wp)                        , INTENT(in) ::   p2dt   ! time step
2063      !
2064      INTEGER  :: ji, jj, jk ! dummy loop indices
2065      REAL(wp) :: zr_tscale
2066      REAL(wp) :: za1, za2, za3, za4, zdiff
2067      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt
2068      !!---------------------------------------------------------------------
2069      !
2070      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_zdf')
2071      !
2072      zr_tscale = 0.5
2073      !
2074      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt) 
2075      !
2076      !
2077
2078      zwt(:,:,1:jpk) = 1.e-10
2079
2080      zwt(:,:,1) = 0.0_wp
2081!      DO jk = 2, jpkm1
2082!         DO jj = 1, jpj
2083!            DO ji = 1, jpi
2084!               zwt(ji,jj,jk) = zwt(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2085!            END DO
2086!         END DO
2087!      END DO
2088      !
2089      !
2090      ! Set diffusivity (homogeneous to an inverse time scale)
2091      !
2092      DO jk = 2, jpkm1
2093         DO jj = 2, jpjm1
2094            DO ji = fs_2, fs_jpim1
2095! Taper a little bit:
2096                za1 = tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)
2097                za2 = tilde_e3t_n(ji,jj,jk  )+e3t_0(ji,jj,jk  )
2098                za4 = 0.5_wp * (e3t_0(ji,jj,jk-1) + e3t_0(ji,jj,jk  ))
2099                za3 = 0.5_wp * (za1 + za2) 
2100                zdiff = ABS(za3-za4)/za4
2101!                IF (zdiff>=0.8) THEN
2102!                   zwt(ji,jj,jk) =  zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk)
2103                   zwt(ji,jj,jk) =  dsm(ji,jj)/ht_0(ji,jj)*(1._wp-tanh((mbkt(ji,jj)+1-jk)*0.2))*tmask(ji,jj,jk)
2104
2105!                ELSE
2106!                   zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk)
2107!                ENDIF
2108            END DO
2109         END DO
2110      END DO
2111      !
2112      !
2113      ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
2114      DO jk = 1, jpkm1
2115         DO jj = 2, jpjm1
2116            DO ji = fs_2, fs_jpim1   ! vector opt.
2117                zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / fse3w(ji,jj,jk  )
2118                zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / fse3w(ji,jj,jk+1)
2119                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
2120            END DO
2121         END DO
2122      END DO
2123      !
2124      !! Matrix inversion from the first level
2125      !!----------------------------------------------------------------------
2126      DO jj = 2, jpjm1
2127         DO ji = fs_2, fs_jpim1
2128            zwt(ji,jj,1) = zwd(ji,jj,1)
2129         END DO
2130      END DO
2131      DO jk = 2, jpkm1
2132         DO jj = 2, jpjm1
2133            DO ji = fs_2, fs_jpim1
2134               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
2135            END DO
2136         END DO
2137      END DO
2138      !         
2139      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
2140      DO jk = 2, jpkm1
2141         DO jj = 2, jpjm1
2142            DO ji = fs_2, fs_jpim1
2143               tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * tilde_e3t_a(ji,jj,jk-1)
2144            END DO
2145         END DO
2146      END DO
2147      ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
2148      DO jj = 2, jpjm1
2149         DO ji = fs_2, fs_jpim1
2150            tilde_e3t_a(ji,jj,jpkm1) = tilde_e3t_a(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
2151         END DO
2152      END DO
2153      DO jk = jpk-2, 1, -1
2154         DO jj = 2, jpjm1
2155            DO ji = fs_2, fs_jpim1
2156               tilde_e3t_a(ji,jj,jk) = ( tilde_e3t_a(ji,jj,jk) - zws(ji,jj,jk) * tilde_e3t_a(ji,jj,jk+1) )   &
2157                  &            / zwt(ji,jj,jk) * tmask(ji,jj,jk)
2158            END DO
2159         END DO
2160      END DO
2161      !
2162      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt) 
2163      !
2164      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_zdf')
2165      !
2166   END SUBROUTINE dom_vvl_zdf
2167
2168   SUBROUTINE dom_vvl_zdf2( kt, p2dt ) 
2169      !!----------------------------------------------------------------------
2170      !!                  ***  ROUTINE dom_vvl_zdf  ***
2171      !!
2172      !! ** Purpose :   Do vertical interface diffusion
2173      !!
2174      !! ** Method  : 
2175      !!
2176      !! ** Action  :
2177      !!---------------------------------------------------------------------
2178      USE oce     , ONLY:   zwd => ua       , zws => va         ! (ua,va) used as 3D workspace
2179      !
2180      INTEGER                         , INTENT(in) ::   kt     ! ocean time-step index
2181      REAL(wp)                        , INTENT(in) ::   p2dt   ! time step
2182      !
2183      INTEGER  :: ji, jj, jk, kbot, kbotm1 ! dummy loop indices
2184      REAL(wp) :: zr_tscale
2185      REAL(wp) :: za1, za2, za3, za4, zdiff
2186      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt, zw
2187      !!---------------------------------------------------------------------
2188      !
2189      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_zdf')
2190      !
2191      zr_tscale = 0.5
2192      !
2193      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zw) 
2194      !
2195      ! Compute internal interfaces depths:
2196      zw(:,:,1) = 0._wp
2197      DO jk = 2, jpk
2198         DO jj = 2, jpjm1
2199            DO ji = fs_2, fs_jpim1   ! vector opt.
2200               zw(ji,jj,jk) = zw(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))  * tmask(ji,jj,jk-1)
2201            END DO
2202         END DO
2203      END DO
2204      !
2205      ! Set diffusivities at interfaces
2206      zwt(:,:,:) = 0.00000001_wp * tmask(:,:,:)     
2207      zwt(:,:,1) = 0._wp
2208      !
2209      !
2210      ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
2211      DO jk = 2, jpkm1
2212         DO jj = 2, jpjm1
2213            DO ji = fs_2, fs_jpim1   ! vector opt.
2214                zwi(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk-1) + zwt(ji,jj,jk  )) & 
2215                              &                   / fse3w(ji,jj,jk-1) / fse3t(ji,jj,jk-1) &
2216                              &                   * ht(ji,jj) 
2217                zws(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk  ) + zwt(ji,jj,jk+1)) & 
2218                              &                   / fse3w(ji,jj,jk  ) / fse3t(ji,jj,jk  ) &
2219                              &                   * ht(ji,jj) 
2220
2221                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
2222            END DO
2223         END DO
2224      END DO
2225      ! Boundary conditions (Neumann)
2226      DO jj = 2, jpjm1
2227         DO ji = fs_2, fs_jpim1
2228            zwi(ji,jj,1) = 0._wp
2229            zws(ji,jj,1) = 0._wp
2230            zwd(ji,jj,1) = 1._wp
2231            zw (ji,jj,1) = 0._wp
2232            !
2233            zwd(ji,jj,2) = zwd(ji,jj,2) +  zwi(ji,jj,2)
2234            zwi(ji,jj,2) = 0._wp
2235!            zwi(ji,jj,2) = 0._wp
2236!            zws(ji,jj,2) = 0._wp
2237!            zwd(ji,jj,2) = 1._wp
2238         END DO
2239      END DO
2240      DO jj = 2, jpjm1
2241         DO ji = fs_2, fs_jpim1
2242            kbot   = mbkt(ji,jj) + 1
2243            kbotm1 = mbkt(ji,jj)
2244            zwi(ji,jj,kbot  ) = 0._wp
2245            zws(ji,jj,kbot  ) = 0._wp
2246            zwd(ji,jj,kbot  ) = 1._wp
2247            !
2248            zwd(ji,jj,kbotm1) = zwd(ji,jj,kbotm1) +  zws(ji,jj,kbotm1)
2249            zws(ji,jj,kbotm1) = 0._wp
2250!            zwi(ji,jj,kbotm1) = 0._wp
2251!            zws(ji,jj,kbotm1) = 0._wp
2252!            zwd(ji,jj,kbotm1) = 1._wp
2253         END DO
2254      END DO
2255      !
2256      DO jk = 2, jpkm1
2257         DO jj = 2, jpjm1
2258            DO ji = fs_2, fs_jpim1
2259               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
2260            END DO
2261         END DO
2262      END DO
2263      !
2264      DO jk = 2, jpk
2265         DO jj = 2, jpjm1
2266            DO ji = fs_2, fs_jpim1    ! vector opt.
2267               zwi(ji,jj,jk) = zw(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) *zwi(ji,jj,jk-1)
2268            END DO
2269         END DO
2270      END DO
2271      !
2272      DO jk = jpk-1, 2, -1
2273         DO jj = 2, jpjm1
2274            DO ji = fs_2, fs_jpim1    ! vector opt.
2275               zw(ji,jj,jk) = ( zwi(ji,jj,jk) - zws(ji,jj,jk) * zw(ji,jj,jk+1) ) / zwd(ji,jj,jk)
2276            END DO
2277         END DO
2278      END DO
2279      !
2280      ! Revert to thicknesses anomalies:
2281      DO jk = 1, jpkm1
2282         DO jj = 2, jpjm1
2283            DO ji = fs_2, fs_jpim1   ! vector opt.
2284               tilde_e3t_a(ji,jj,jk) = (zw(ji,jj,jk+1)-zw(ji,jj,jk)-e3t_0(ji,jj,jk))* tmask(ji,jj,jk)
2285            END DO
2286         END DO
2287      END DO
2288      !
2289      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zw) 
2290      !
2291      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_zdf')
2292      !
2293   END SUBROUTINE dom_vvl_zdf2
2294
2295   SUBROUTINE dom_vvl_regrid( kt )
2296      !!----------------------------------------------------------------------
2297      !!                ***  ROUTINE dom_vvl_regrid  ***
2298      !!                   
2299      !! ** Purpose :  Ensure "well-behaved" vertical grid
2300      !!
2301      !! ** Method  :  More or less adapted from references below.
2302      !!
2303      !! ** Action  :  Ensure that thickness are above a given value, spaced enough
2304      !!               and revert to Eulerian coordinates near the bottom.     
2305      !!
2306      !! References :  Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction
2307      !!               with a Model Combining Terrain-following and Isentropic
2308      !!               coordinates. Part I: Model Description. Monthly Weather Rev.,
2309      !!               121, 1770-1785.
2310      !!               Toy, M., 2011: Incorporating Condensational Heating into a
2311      !!               Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic-
2312      !!               Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954.
2313      !!----------------------------------------------------------------------
2314      !! * Arguments
2315      INTEGER, INTENT( in )               :: kt         ! time step
2316
2317      !! * Local declarations
2318      INTEGER                             :: ji, jj, jk ! dummy loop indices
2319      LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond
2320      LOGICAL                             :: ll_e3tsurf_const, ll_zdiff_cond, ll_blpdiff_cond
2321      INTEGER                             :: jkbot
2322      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd, dzmin_vvl
2323      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj
2324      REAL(wp)                            :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn
2325      REAL(wp)                            :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot
2326      REAL(wp)                            :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim
2327      REAL(wp), POINTER, DIMENSION(:,:)   :: zdw
2328      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zwdw_b
2329      !!----------------------------------------------------------------------
2330
2331      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_regrid')
2332      !
2333      CALL wrk_alloc( jpi, jpj, jpk, zwdw) 
2334      !
2335      ! Some user defined parameters below:
2336      ll_chk_bot2top   = .TRUE.
2337      ll_chk_top2bot   = .TRUE.
2338      ll_e3tsurf_const = .FALSE.
2339      dzmin_vvl  = 1._wp   ! Absolute minimum depth (in meters)
2340      zfrch_stp  = 0.5_wp   ! Maximum fractionnal thickness change in one time step (<= 1.)
2341      zfrch_rel  = 0.0_wp   ! Maximum relative thickness change in the vertical (<= 1.)
2342      zfrac_bot  = 0.01_wp  ! Fraction of bottom level allowed to change
2343      zscal_bot  = 2._wp    ! Depth lengthscale
2344      ll_zdiff_cond    = .FALSE.  ! Conditionnal vertical diffusion of interfaces
2345         zvdiff  = 0.1_wp   ! m/s
2346         zvlim   = 0.4_wp   ! max d2h/dh
2347      ll_lapdiff_cond  = .FALSE.  ! Conditionnal Laplacian diffusion of interfaces
2348         zhdiff  = 0.1_wp   ! ad.
2349         zhlim   = 0.01_wp  ! max lap(z)/e1
2350      ll_blpdiff_cond  = .FALSE.  ! Conditionnal Bilaplacian diffusion of interfaces
2351         zhdiff2 = 0.8_wp   ! ad.
2352         zhlim2  = 5.e-11_wp  ! max bilap(z)/e1**3
2353      ! ---------------------------------------------------------------------------------------
2354      !
2355      ! Set arrays determining maximum vertical displacement at the bottom:
2356      !--------------------------------------------------------------------
2357      IF ( kt==nit000 ) THEN
2358         DO jj = 2, jpjm1
2359            DO ji = 2, jpim1
2360               jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1))
2361               i_int_bot(ji,jj) = jk
2362            END DO
2363         END DO
2364         dsm(:,:) = REAL( i_int_bot(:,:), wp )   ;   CALL lbc_lnk(dsm(:,:),'T',1.)
2365         i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 )
2366
2367         CALL wrk_alloc( jpi, jpj, zdw )
2368         DO jj = 2, jpjm1
2369            DO ji = 2, jpim1
2370               zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji  ,jj,1), &
2371                          &     ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), &
2372                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj  ,1), &
2373                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1)  )
2374               zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall ) 
2375            END DO
2376         END DO
2377         CALL lbc_lnk( zdw(:,:), 'T', 1. )
2378
2379         DO jj = 2, jpjm1
2380            DO ji = 2, jpim1
2381               dsm(ji,jj) = 1._wp/16._wp * (            zdw(ji-1,jj-1) + zdw(ji+1,jj-1)      &
2382                  &                           +         zdw(ji-1,jj+1) + zdw(ji+1,jj+1)      &
2383                  &                           + 2._wp*( zdw(ji  ,jj-1) + zdw(ji-1,jj  )      &
2384                  &                           +         zdw(ji+1,jj  ) + zdw(ji  ,jj+1) )    &
2385                  &                           + 4._wp*  zdw(ji  ,jj  )                       )
2386            END DO
2387         END DO
2388
2389         CALL lbc_lnk( dsm(:,:), 'T', 1. )
2390         CALL wrk_dealloc( jpi, jpj, zdw)         
2391     
2392         IF (ln_zps) THEN
2393            DO jj = 1, jpj
2394               DO ji = 1, jpi
2395                  jk = i_int_bot(ji,jj)
2396                  hsm(ji,jj) = zfrac_bot * e3w_1d(jk)
2397               END DO
2398            END DO
2399         ELSE
2400            DO jj = 1, jpj
2401               DO ji = 1, jpi
2402                  jk = i_int_bot(ji,jj)
2403                  hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk)
2404               END DO
2405            END DO
2406         ENDIF
2407      END IF
2408
2409      ! Provisionnal interface depths:
2410      !-------------------------------
2411      zwdw(:,:,1) = 0.e0
2412      DO jj = 1, jpj
2413         DO ji = 1, jpi
2414            DO jk = 2, jpk
2415               zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + & 
2416                              & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2417            END DO
2418         END DO
2419      END DO
2420      !
2421      ! Conditionnal horizontal Laplacian diffusion:
2422      !---------------------------------------------
2423      IF ( ll_lapdiff_cond ) THEN
2424         CALL wrk_alloc( jpi, jpj, jpk, zwdw_b)
2425         !
2426         zwdw_b(:,:,1) = 0._wp
2427         DO jj = 1, jpj
2428            DO ji = 1, jpi
2429               DO jk=2,jpk
2430                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
2431                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2432               END DO
2433            END DO
2434         END DO
2435         !
2436         DO jk = 2, jpkm1
2437            DO jj = 1, jpjm1
2438               DO ji = 1, fs_jpim1   ! vector opt.                 
2439                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
2440                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
2441                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
2442                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
2443               END DO
2444            END DO
2445         END DO
2446           
2447         DO jk = 2, jpkm1
2448            DO jj = 2, jpjm1
2449               DO ji = fs_2, fs_jpim1   ! vector opt.
2450                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
2451                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
2452                  zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e12t(ji,jj)), 0._wp)
2453                  ztmp = SIGN(zh2, ztmp1)
2454                  zeu2 = zhdiff * e12t(ji,jj)*e12t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj))
2455                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
2456               END DO
2457            END DO
2458         END DO         
2459         !
2460         CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b)
2461         !
2462      ENDIF
2463
2464      ! Conditionnal horizontal Bilaplacian diffusion:
2465      !-----------------------------------------------
2466      IF ( ll_blpdiff_cond ) THEN
2467         CALL wrk_alloc( jpi, jpj, jpk, zwdw_b)
2468         !
2469         zwdw_b(:,:,1) = 0._wp
2470         DO jj = 1, jpj
2471            DO ji = 1, jpi
2472               DO jk=2,jpk
2473                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
2474                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2475               END DO
2476            END DO
2477         END DO
2478         !
2479         DO jk = 2, jpkm1
2480            DO jj = 1, jpjm1
2481               DO ji = 1, fs_jpim1   ! vector opt.                 
2482                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
2483                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
2484                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
2485                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
2486               END DO
2487            END DO
2488         END DO
2489           
2490         DO jk = 2, jpkm1
2491            DO jj = 2, jpjm1
2492               DO ji = fs_2, fs_jpim1   ! vector opt.
2493                  zwdw_b(ji,jj,jk) = -( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
2494                                &  +    (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
2495               END DO
2496            END DO
2497         END DO   
2498         !
2499         CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )     
2500         !
2501         DO jk = 2, jpkm1
2502            DO jj = 1, jpjm1
2503               DO ji = 1, fs_jpim1   ! vector opt.                 
2504                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) &
2505                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
2506                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
2507                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
2508               END DO
2509            END DO
2510         END DO
2511         !   
2512         DO jk = 2, jpkm1
2513            DO jj = 2, jpjm1
2514               DO ji = fs_2, fs_jpim1   ! vector opt.
2515                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    &
2516                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
2517                  zh2 = MAX(abs(ztmp1)-zhlim2, 0._wp)
2518                  ztmp = SIGN(zh2, ztmp1)
2519                  zeu2 = zhdiff2 * e12t(ji,jj)*e12t(ji,jj) / 16._wp
2520                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
2521               END DO
2522            END DO
2523         END DO   
2524         !
2525         CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b)
2526         !
2527      ENDIF
2528
2529      ! Conditionnal vertical diffusion:
2530      !---------------------------------
2531      IF ( ll_zdiff_cond ) THEN
2532         DO jk = 2, jpkm1
2533            DO jj = 2, jpjm1
2534               DO ji = fs_2, fs_jpim1   ! vector opt.   
2535                  ztmp  = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) & 
2536                            -(tilde_e3t_b(ji,jj,jk  )+e3t_0(ji,jj,jk  ))*tmask(ji,jj,jk  ) ) 
2537                  ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) & 
2538                        &           +tilde_e3t_b(ji,jj,jk  ) + e3t_0(ji,jj,jk  ) )
2539                  zh2   = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp)
2540                  ztmp  = SIGN(zh2, ztmp)
2541                  IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0
2542                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk)
2543               END DO
2544            END DO
2545         END DO
2546      ENDIF
2547      !
2548      ! Check grid from the bottom to the surface
2549      !------------------------------------------
2550      IF ( ll_chk_bot2top ) THEN
2551         DO jj = 2, jpjm1
2552            DO ji = 2, jpim1
2553               jkbot = mbkt(ji,jj)                   
2554               DO jk = jkbot,2,-1
2555                  !
2556                  zh_0   = e3t_0(ji,jj,jk)
2557                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1))
2558                  zh_old = tilde_e3t_a(ji,jj,jk  ) + zh_0
2559!                  zh_dwn = tilde_e3t_a(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)
2560                  zh_min = MIN(zh_0/3._wp, dzmin_vvl)
2561                  !
2562                  ! Set maximum and minimum vertical excursions   
2563                  ztmph = hsm(ji,jj)
2564                  ztmpd = dsm(ji,jj)
2565                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)               
2566                  zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 )
2567                  zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff)
2568                  zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 )
2569                  zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff)
2570                  !
2571                  ! New layer thickness:
2572                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
2573                  !
2574                  ! Ensure minimum layer thickness:                 
2575!                  zh_new = MIN(zh_old, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) )
2576                  zh_new = cush(zh_new, zh_min)
2577                  !
2578                  ! Final flux:
2579                  zdiff = (zh_new - zh_old) * tmask(ji,jj,jk)
2580                  !
2581                  ! Limit flux:
2582                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
2583                  zdiff = SIGN(ztmp, zh_new - zh_old)
2584                  zh_new = zdiff + zh_old
2585                  !
2586                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
2587                  zwdw(ji,jj,jk)          = zwdw(ji,jj,jk+1) - zh_new
2588                  tilde_e3t_a(ji,jj,jk-1) = (-zdiff + tilde_e3t_a(ji,jj,jk-1) ) * tmask(ji,jj,jk-1)           
2589               END DO   
2590            END DO
2591         END DO
2592      END IF   
2593      !
2594      ! Check grid from the surface to the bottom
2595      !------------------------------------------
2596      IF ( ll_chk_top2bot ) THEN     
2597         DO jj = 2, jpjm1
2598            DO ji = 2, jpim1
2599               jkbot = mbkt(ji,jj)   
2600               DO jk = 1, jkbot-1
2601                  !
2602                  zh_0   = e3t_0(ji,jj,jk)
2603                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1))
2604                  zh_old = tilde_e3t_a(ji,jj,jk  ) + zh_0   
2605!                  zh_up  = tilde_e3t_a(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)
2606                  zh_min = MIN(zh_0/3._wp, dzmin_vvl)
2607                  IF ((jk<=5).AND.ll_e3tsurf_const) zh_min = MAX(e3t_0(ji,jj,1)/3._wp, dzmin_vvl)
2608                  !
2609                  ! Ensure minimum layer thickness:
2610!                  zh_new=MIN(zh_old, zh_up * zfrch_rel / (2._wp-zfrch_rel) )
2611                  zh_new = cush(zh_old, zh_min)
2612                  !
2613                  ! Final flux:
2614                  zdiff = (zh_new -zh_old) * tmask(ji,jj,jk)
2615                  !
2616                  ! Limit flux:                 
2617                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
2618                  zdiff = SIGN(ztmp, zdiff)
2619                  zh_new = zdiff + zh_old
2620                  !
2621                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
2622                  zwdw(ji,jj,jk+1)        = zwdw(ji,jj,jk) + zh_new
2623                  tilde_e3t_a(ji,jj,jk+1) = (-zdiff + tilde_e3t_a(ji,jj,jk+1) ) * tmask(ji,jj,jk+1)
2624               END DO
2625               !               
2626            END DO
2627         END DO
2628      ENDIF
2629      !
2630      CALL wrk_dealloc( jpi, jpj, jpk, zwdw ) 
2631      !
2632      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_regrid')
2633      !
2634   END SUBROUTINE dom_vvl_regrid
2635
2636   FUNCTION cush(hin, hmin)  RESULT(hout)
2637      !!----------------------------------------------------------------------
2638      !!                 ***  FUNCTION cush  ***
2639      !!
2640      !! ** Purpose :   
2641      !!
2642      !! ** Method  :
2643      !!
2644      !!----------------------------------------------------------------------
2645      IMPLICIT NONE
2646      REAL(wp), INTENT(in) ::  hin, hmin
2647      REAL(wp)             ::  hout, zx, zh_cri
2648      !!----------------------------------------------------------------------
2649      zh_cri = 3._wp * hmin
2650      !
2651      IF ( hin<=0._wp ) THEN
2652         hout = hmin
2653      !
2654      ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN
2655         zx = hin/zh_cri
2656         hout = hmin * (1._wp + zx + zx*zx)     
2657      !
2658      ELSEIF ( hin>zh_cri ) THEN
2659         hout = hin
2660      !
2661      ENDIF
2662      !
2663   END FUNCTION cush
2664
2665   FUNCTION cush_max(hin, hmax)  RESULT(hout)
2666      !!----------------------------------------------------------------------
2667      !!                 ***  FUNCTION cush  ***
2668      !!
2669      !! ** Purpose :   
2670      !!
2671      !! ** Method  :
2672      !!
2673      !!----------------------------------------------------------------------
2674      IMPLICIT NONE
2675      REAL(wp), INTENT(in) ::  hin, hmax
2676      REAL(wp)             ::  hout, hmin, zx, zh_cri
2677      !!----------------------------------------------------------------------
2678      hmin = 0.1_wp * hmax
2679      zh_cri = 3._wp * hmin
2680      !
2681      IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN
2682         zx = (hmax-hin)/zh_cri
2683         hout = hmax - hmin * (1._wp + zx + zx*zx)
2684         !
2685      ELSEIF ( hin>(hmax-zh_cri) ) THEN
2686         hout = hmax - hmin
2687         !
2688      ELSE
2689         hout = hin
2690         !
2691      ENDIF
2692      !
2693   END FUNCTION cush_max
2694
2695   SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin )
2696      !!----------------------------------------------------------------------
2697      !!                  ***  ROUTINE dom_vvl_adv_fct  ***
2698      !!
2699      !! **  Purpose :  Do thickness advection
2700      !!
2701      !! **  Method  :  FCT scheme to ensure positivity
2702      !!
2703      !! **  Action  : - Update pta thickness tendency and diffusive fluxes
2704      !!               - this is the total trend, hence it does include sea level motions
2705      !!               - Upstream corrections to antidiffusive fluxes ensure
2706      !!                 that barotropic transport matches what is contained in input fluxes
2707      !!----------------------------------------------------------------------
2708      !
2709      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
2710      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta    ! thickness baroclinic trend
2711      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)    ::   uin, vin  ! input velocities
2712      !
2713      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
2714      INTEGER  ::   ikbu, ikbv, ibot
2715      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
2716      REAL(wp) ::   zdi, zdj, zmin           !   -      -
2717      REAL(wp) ::   zfp_ui, zfp_vj           !   -      -
2718      REAL(wp) ::   zfm_ui, zfm_vj           !   -      -
2719      REAL(wp) ::   zfp_hi, zfp_hj           !   -      -
2720      REAL(wp) ::   zfm_hi, zfm_hj           !   -      -
2721      REAL(wp), POINTER, DIMENSION(:,:)   :: zbu, zbv, zhu_b, zhv_b
2722      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwi
2723      !!----------------------------------------------------------------------
2724      !
2725      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_adv_fct')
2726      !
2727      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
2728      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwi)
2729      !
2730      ! 1. Initializations
2731      ! ------------------
2732      !
2733      IF( neuler == 0 .AND. kt == nit000 ) THEN
2734         z2dtt =  rdt
2735      ELSE
2736         z2dtt = 2.0_wp * rdt
2737      ENDIF
2738      !
2739      zwi(:,:,:) = 0.e0
2740      zwx(:,:,:) = 0.e0 
2741      zwy(:,:,:) = 0.e0
2742      !
2743      !
2744      ! 2. upstream advection with initial mass fluxes & intermediate update
2745      ! --------------------------------------------------------------------
2746      DO jk = 1, jpkm1
2747         DO jj = 1, jpjm1
2748            DO ji = 1, fs_jpim1   ! vector opt.
2749               !
2750               zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
2751               zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
2752               zfp_hi = fse3t_b(ji  ,jj  ,jk)
2753               zfm_hi = fse3t_b(ji+1,jj  ,jk)
2754               zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2755               !
2756               zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
2757               zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
2758               zfp_hj = fse3t_b(ji  ,jj  ,jk)
2759               zfm_hj = fse3t_b(ji  ,jj+1,jk)
2760               zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2761            END DO
2762         END DO
2763      END DO
2764
2765      IF ( .NOT.ln_sco ) THEN
2766         ! Correct bottom upstream fluxes
2767         ! considering "shelf horizon depths"
2768         DO jj = 1, jpjm1
2769            DO ji = 1, fs_jpim1   ! vector opt.
2770               ! upstream scheme
2771               ikbu = mbku(ji,jj)
2772               ikbv = mbkv(ji,jj)
2773               zfp_ui = uin(ji,jj,ikbu) + ABS( uin(ji,jj,ikbu) )
2774               zfm_ui = uin(ji,jj,ikbu) - ABS( uin(ji,jj,ikbu) )
2775               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
2776               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
2777               zwx(ji,jj,ikbu) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,ikbu)
2778               !
2779               zfp_vj = vin(ji,jj,ikbv) + ABS( vin(ji,jj,ikbv) )
2780               zfm_vj = vin(ji,jj,ikbv) - ABS( vin(ji,jj,ikbv) )
2781               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
2782               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
2783               zwy(ji,jj,ikbv) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,ikbv)
2784            END DO
2785         END DO
2786      ENDIF
2787
2788      ! total advective trend
2789      DO jk = 1, jpkm1
2790         DO jj = 2, jpjm1
2791            DO ji = fs_2, fs_jpim1   ! vector opt.
2792               zbtr = r1_e12t(ji,jj)
2793               ! total intermediate advective trends
2794               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )  &
2795                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )  )
2796               !
2797               ! update and guess with monotonic sheme
2798               pta(ji,jj,jk) =   pta(ji,jj,jk) + ztra
2799               zwi(ji,jj,jk) = (fse3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk)
2800            END DO
2801         END DO
2802      END DO
2803      !                             ! Lateral boundary conditions on zwi  (unchanged sign)
2804      CALL lbc_lnk( zwi, 'T', 1. ) 
2805
2806      IF ( ln_vvl_dbg ) THEN
2807         zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 ) 
2808         IF( lk_mpp )   CALL mpp_min( zmin )
2809         IF( zmin < 0._wp) THEN
2810            IF(lwp) CALL ctl_stop('vvl_adv: CFL issue here')
2811         ENDIF
2812      ENDIF
2813
2814      ! 3. antidiffusive flux : high order minus low order
2815      ! --------------------------------------------------
2816      ! antidiffusive flux on i and j
2817      DO jk = 1, jpkm1
2818         DO jj = 1, jpjm1
2819            DO ji = 1, fs_jpim1   ! vector opt.
2820               zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * fse3u_n(ji,jj,jk) &
2821                                & - zwx(ji,jj,jk)) * umask(ji,jj,jk)
2822               zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * fse3v_n(ji,jj,jk) & 
2823                                & - zwy(ji,jj,jk)) * vmask(ji,jj,jk)
2824               !
2825               ! Update advective fluxes
2826               un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk)
2827               vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk)
2828            END DO
2829         END DO
2830      END DO
2831     
2832      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
2833
2834      ! 4. monotonicity algorithm
2835      ! -------------------------
2836      CALL nonosc_2d( fse3t_b(:,:,:), zwx, zwy, zwi, z2dtt )
2837
2838      ! 5. final trend with corrected fluxes
2839      ! ------------------------------------
2840      !
2841      ! Update advective fluxes
2842      un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:)
2843      vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:)
2844      !
2845      DO jk = 1, jpkm1
2846         DO jj = 2, jpjm1
2847            DO ji = fs_2, fs_jpim1   ! vector opt. 
2848               !
2849               zbtr = r1_e12t(ji,jj)
2850               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
2851                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
2852               ! add them to the general tracer trends
2853               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
2854            END DO
2855         END DO
2856      END DO
2857
2858      !
2859      ! 6. Correct barotropic flux
2860      ! --------------------------
2861      ! Compute barotropic flux difference:
2862      zbu(:,:) = 0.e0
2863      zbv(:,:) = 0.e0
2864      DO jj = 1, jpj
2865         DO ji = 1, jpi   ! vector opt.
2866            DO jk = 1, jpkm1
2867               zbu(ji,jj) = zbu(ji,jj) - un_td(ji,jj,jk) * umask(ji,jj,jk)
2868               zbv(ji,jj) = zbv(ji,jj) - vn_td(ji,jj,jk) * vmask(ji,jj,jk)
2869            END DO
2870         END DO
2871      ENDDO
2872
2873      ! Compute upstream depths:
2874      zhu_b(:,:) = 0.e0
2875      zhv_b(:,:) = 0.e0
2876      IF ( ln_sco ) THEN; ibot=0 ; ELSE ; ibot=1 ; ENDIF
2877
2878      DO jj = 1, jpjm1
2879         DO ji = 1, jpim1   ! vector opt.
2880            zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
2881            ikbu = mbku(ji,jj) - ibot 
2882            DO jk = 1, ikbu
2883               zfp_hi = fse3t_b(ji  ,jj  ,jk)
2884               zfm_hi = fse3t_b(ji+1,jj  ,jk)
2885               zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi &
2886                            &                  + (1._wp-zdi) * zfm_hi &
2887                            &                ) * umask(ji,jj,jk)
2888            END DO
2889            !
2890            zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
2891            ikbv = mbkv(ji,jj) - ibot
2892            DO jk = 1, ikbv
2893               zfp_hj = fse3t_b(ji  ,jj  ,jk)
2894               zfm_hj = fse3t_b(ji  ,jj+1,jk)
2895               zhv_b(ji,jj) = zhv_b(ji,jj) + (          zdj  * zfp_hj &
2896                            &                  + (1._wp-zdj) * zfm_hj &
2897                            &                ) * vmask(ji,jj,jk)
2898            END DO
2899         END DO
2900      END DO
2901
2902      IF ( .NOT.ln_sco ) THEN
2903         ! Correct bottom value
2904         ! considering "shelf horizon depth"
2905         DO jj = 1, jpjm1
2906            DO ji = 1, jpim1   ! vector opt.
2907               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
2908               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
2909               ikbu = mbku(ji,jj)
2910               ikbv = mbkv(ji,jj)
2911               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
2912               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
2913               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
2914               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
2915               zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi &
2916                            &                 + (1._wp-zdi) * zfm_hi &
2917                            &                ) * umask(ji,jj,ikbu)
2918               zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj &
2919                            &                 + (1._wp-zdj) * zfm_hj &
2920                            &                ) * vmask(ji,jj,ikbv)
2921            END DO
2922         END DO
2923      ENDIF
2924
2925      ! Corrective barotropic velocity (times hor. scale factor)
2926      zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1))
2927      zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1))
2928
2929      CALL lbc_lnk( zbu(:,:), 'U', -1. )
2930      CALL lbc_lnk( zbv(:,:), 'V', -1. )
2931      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
2932     
2933      ! Set corrective fluxes in upstream direction:
2934      !
2935      zwx(:,:,:) = 0.e0
2936      zwy(:,:,:) = 0.e0
2937      DO jj = 1, jpjm1
2938         DO ji = 1, fs_jpim1   ! vector opt.
2939            ! upstream scheme
2940            zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
2941            zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
2942            zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
2943            zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
2944            DO jk = 1, jpkm1
2945               zfp_hi = fse3t_b(ji  ,jj  ,jk)
2946               zfm_hi = fse3t_b(ji+1,jj  ,jk)
2947               !
2948               zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2949
2950               zfp_hj = fse3t_b(ji  ,jj  ,jk)
2951               zfm_hj = fse3t_b(ji  ,jj+1,jk)
2952               !
2953               zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2954            END DO
2955         END DO
2956      END DO
2957
2958      IF ( .NOT.ln_sco ) THEN
2959      ! Bottom correction:
2960      DO jj = 1, jpjm1
2961         DO ji = 1, fs_jpim1   ! vector opt.
2962            ! upstream scheme
2963            ikbu = mbku(ji,jj)
2964            ikbv = mbkv(ji,jj)
2965            !
2966            zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
2967            zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
2968            zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
2969            zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
2970            !
2971            zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp)
2972            zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp)
2973            !
2974            zwx(ji,jj,ikbu) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi )
2975            !
2976            zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp)
2977            zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp)
2978            !
2979            zwy(ji,jj,ikbv) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj )
2980         END DO
2981      END DO
2982      ENDIF
2983
2984      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions
2985
2986      un_td(:,:,:) = un_td(:,:,:) + zwx(:,:,:)
2987      vn_td(:,:,:) = vn_td(:,:,:) + zwy(:,:,:)
2988      !
2989      ! Update trend with corrective fluxes:
2990      DO jk = 1, jpkm1
2991         DO jj = 2, jpjm1
2992            DO ji = fs_2, fs_jpim1   ! vector opt. 
2993               !
2994               zbtr = r1_e12t(ji,jj)
2995               ! total advective trends
2996               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
2997                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
2998               ! add them to the general tracer trends
2999               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
3000            END DO
3001         END DO
3002      END DO
3003
3004      CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
3005      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi)
3006      !
3007      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_adv_fct')
3008      !
3009   END SUBROUTINE dom_vvl_adv_fct
3010
3011   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt )
3012      !!---------------------------------------------------------------------
3013      !!                    ***  ROUTINE nonosc_2d  ***
3014      !!     
3015      !! **  Purpose :   compute monotonic thickness fluxes from the upstream
3016      !!       scheme and the before field by a nonoscillatory algorithm
3017      !!
3018      !! **  Method  :   ... ???
3019      !!       warning : pbef and paft must be masked, but the boundaries
3020      !!       conditions on the fluxes are not necessary zalezak (1979)
3021      !!       drange (1995) multi-dimensional forward-in-time and upstream-
3022      !!       in-space based differencing for fluid
3023      !!----------------------------------------------------------------------
3024      !
3025      !!----------------------------------------------------------------------
3026      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
3027      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
3028      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb        ! monotonic fluxes in the 3 directions
3029      !
3030      INTEGER ::   ji, jj, jk   ! dummy loop indices
3031      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
3032      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
3033      REAL(wp) ::   zupip1, zupim1, zupjp1, zupjm1, zupb, zupa
3034      REAL(wp) ::   zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa
3035      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
3036      !!----------------------------------------------------------------------
3037      !
3038      IF( nn_timing == 1 )  CALL timing_start('nonosc2')
3039      !
3040      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
3041      !
3042
3043      zbig  = 1.e+40_wp
3044      zrtrn = 1.e-15_wp
3045      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp
3046
3047
3048      ! Search local extrema
3049      ! --------------------
3050      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
3051      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   &
3052         &        paft * tmask - zbig * ( 1.e0 - tmask )  )
3053      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   &
3054         &        paft * tmask + zbig * ( 1.e0 - tmask )  )
3055
3056      DO jk = 1, jpkm1
3057         z2dtt = p2dt
3058         DO jj = 2, jpjm1
3059            DO ji = fs_2, fs_jpim1   ! vector opt.
3060
3061               ! search maximum in neighbourhood
3062               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
3063                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
3064                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ))
3065
3066               ! search minimum in neighbourhood
3067               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
3068                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
3069                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ))
3070
3071               ! positive part of the flux
3072               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
3073                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   
3074
3075               ! negative part of the flux
3076               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
3077                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )
3078
3079               ! up & down beta terms
3080               zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt
3081               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
3082               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
3083            END DO
3084         END DO
3085      END DO
3086
3087      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
3088
3089      ! 3. monotonic flux in the i & j direction (paa & pbb)
3090      ! ----------------------------------------
3091      DO jk = 1, jpkm1
3092         DO jj = 2, jpjm1
3093            DO ji = fs_2, fs_jpim1   ! vector opt.
3094               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
3095               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
3096               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
3097               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
3098
3099               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
3100               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
3101               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
3102               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
3103            END DO
3104         END DO
3105      END DO
3106      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
3107      !
3108      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
3109      !
3110      IF( nn_timing == 1 )  CALL timing_stop('nonosc2')
3111      !
3112   END SUBROUTINE nonosc_2d
3113
3114   !!======================================================================
3115END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.