source: NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90 @ 10164

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

ztilde update, #2126: Add higher order filtering + alternative scale factor decomposition

  • Property svn:keywords set to Id
File size: 120.8 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: vvl option includes z_star and z_tilde coordinates
9   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
10   !!            4.0  !  2018-09  (J. Chanut) improve z_tilde robustness
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
15   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
16   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
17   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
18   !!   dom_vvl_rst      : read/write restart file
19   !!   dom_vvl_ctl      : Check the vvl options
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers
22   USE phycst          ! physical constant
23   USE dom_oce         ! ocean space and time domain
24   USE sbc_oce         ! ocean surface boundary condition
25   USE wet_dry         ! wetting and drying
26   USE usrdef_istate   ! user defined initial state (wad only)
27   USE restart         ! ocean restart
28   !
29   USE in_out_manager  ! I/O manager
30   USE iom             ! I/O manager library
31   USE lib_mpp         ! distributed memory computing library
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33   USE timing          ! Timing
34   USE bdy_oce         ! ocean open boundary conditions
35   USE sbcrnf          ! river runoff
36   USE dynspg_ts, ONLY: un_adv, vn_adv
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC  dom_vvl_init       ! called by domain.F90
42   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
43   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
44   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
45
46   !                                                      !!* Namelist nam_vvl
47   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
49   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
50   LOGICAL          :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
51   LOGICAL          :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
52   LOGICAL          :: ln_vvl_zstar_on_shelf  = .FALSE.    ! revert to zstar on shelves
53   LOGICAL          :: ln_vvl_adv_fct         = .FALSE.    ! Centred thickness advection
54   LOGICAL          :: ln_vvl_adv_cn2         = .TRUE.     ! FCT thickness advection
55   LOGICAL          :: ln_vvl_dbg             = .FALSE.    ! debug control prints
56   LOGICAL          :: ln_vvl_ramp            = .FALSE.    ! Ramp on interfaces displacement
57   LOGICAL          :: ln_vvl_lap             = .FALSE.    ! Laplacian thickness diffusion
58   LOGICAL          :: ln_vvl_blp             = .FALSE.    ! Bilaplacian thickness diffusion
59   LOGICAL          :: ln_vvl_regrid          = .FALSE.    ! ensure layer separation
60   LOGICAL          :: ll_shorizd             = .FALSE.    ! Use "shelf horizon depths"
61   LOGICAL          :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
62   !                                                       ! conservation: not used yet
63   INTEGER          :: nn_filt_order=1
64   REAL(wp)         :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian)
65   REAL(wp)         :: rn_ahe3_blp               ! thickness diffusion coefficient (Bilaplacian)
66   REAL(wp)         :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
67   REAL(wp)         :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
68   REAL(wp)         :: rn_day_ramp               ! Duration of linear ramp  [days]
69   REAL(wp)         :: hsmall=0.01_wp            ! small thickness [m]
70
71   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
72   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: un_lf, vn_lf, hdivn_lf    ! low frequency fluxes and divergence
73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors
74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a                 ! baroclinic scale factors
75   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: tildemask                   ! mask tilde tendency
76   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! restoring period for scale factors
77   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! restoring period for low freq. divergence
78   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: hsm, dsm                    !
79   INTEGER         , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: i_int_bot
80
81   !! * Substitutions
82#  include "vectopt_loop_substitute.h90"
83   !!----------------------------------------------------------------------
84   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
85   !! $Id$
86   !! Software governed by the CeCILL license (see ./LICENSE)
87   !!----------------------------------------------------------------------
88CONTAINS
89
90   INTEGER FUNCTION dom_vvl_alloc()
91      !!----------------------------------------------------------------------
92      !!                ***  FUNCTION dom_vvl_alloc  ***
93      !!----------------------------------------------------------------------
94      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
95      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
96         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
97            &      un_td  (jpi,jpj,jpk)      , vn_td  (jpi,jpj,jpk)     ,                              &
98            &      tildemask(jpi,jpj) , hsm(jpi,jpj) , dsm(jpi,jpj) , i_int_bot(jpi,jpj), STAT = dom_vvl_alloc )
99         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
100         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
101         un_td = 0._wp
102         vn_td = 0._wp
103      ENDIF
104      IF( ln_vvl_ztilde ) THEN
105         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdivn_lf(jpi,jpj,jpk,nn_filt_order),  &
106            &      un_lf(jpi,jpj,jpk,nn_filt_order), vn_lf(jpi,jpj,jpk,nn_filt_order), STAT= dom_vvl_alloc )
107         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
108         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
109      ENDIF
110      !
111   END FUNCTION dom_vvl_alloc
112
113
114   SUBROUTINE dom_vvl_init
115      !!----------------------------------------------------------------------
116      !!                ***  ROUTINE dom_vvl_init  ***
117      !!                   
118      !! ** Purpose :  Initialization of all scale factors, depths
119      !!               and water column heights
120      !!
121      !! ** Method  :  - use restart file and/or initialize
122      !!               - interpolate scale factors
123      !!
124      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
125      !!              - Regrid: e3(u/v)_n
126      !!                        e3(u/v)_b       
127      !!                        e3w_n           
128      !!                        e3(u/v)w_b     
129      !!                        e3(u/v)w_n     
130      !!                        gdept_n, gdepw_n and gde3w_n
131      !!              - h(t/u/v)_0
132      !!              - frq_rst_e3t and frq_rst_hdv
133      !!
134      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
135      !!----------------------------------------------------------------------
136      INTEGER ::   ji, jj, jk
137      INTEGER ::   ii0, ii1, ij0, ij1
138      REAL(wp)::   zcoef, zwgt, ztmp, zhmin, zhmax
139      !!----------------------------------------------------------------------
140      !
141      IF(lwp) WRITE(numout,*)
142      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
143      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
144      !
145      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
146      !
147      !                    ! Allocate module arrays
148      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
149      !
150      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdivn_lf
151      CALL dom_vvl_rst( nit000, 'READ' )
152      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
153      !
154      !                    !== Set of all other vertical scale factors  ==!  (now and before)
155      !                                ! Horizontal interpolation of e3t
156      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
157      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
158      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V
159      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
160      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
161      !                                ! Vertical interpolation of e3t,u,v
162      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
163      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
164      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
165      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
166      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
167      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
168
169      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
170      e3t_a(:,:,:) = e3t_n(:,:,:)
171      e3u_a(:,:,:) = e3u_n(:,:,:)
172      e3v_a(:,:,:) = e3v_n(:,:,:)
173      !
174      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
175      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
176      gdepw_n(:,:,1) = 0.0_wp
177      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
178      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
179      gdepw_b(:,:,1) = 0.0_wp
180      DO jk = 2, jpk                               ! vertical sum
181         DO jj = 1,jpj
182            DO ji = 1,jpi
183               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
184               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
185               !                             ! 0.5 where jk = mikt     
186!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
187               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
188               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
189               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
190                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
191               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
192               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
193               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
194                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
195            END DO
196         END DO
197      END DO
198      !
199      !                    !==  thickness of the water column  !!   (ocean portion only)
200      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
201      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
202      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
203      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
204      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
205      DO jk = 2, jpkm1
206         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
207         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
208         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
209         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
210         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
211      END DO
212      !
213      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
214      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
215      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
216      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
217      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
218
219      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
220      tildemask(:,:) = 1._wp
221
222      IF( ln_vvl_ztilde ) THEN
223!!gm : idea: add here a READ in a file of custumized restoring frequency
224         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
225         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
226         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
227         !
228         IF( ln_vvl_ztilde_as_zstar ) THEN
229            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
230            frq_rst_e3t(:,:) = 0.0_wp 
231            frq_rst_hdv(:,:) = 1.0_wp / rdt
232            tildemask(:,:) = 0._wp
233         ENDIF
234       
235         IF ( ln_vvl_zstar_at_eqtor ) THEN
236            DO jj = 1, jpj
237               DO ji = 1, jpi
238!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
239                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
240                     ! values outside the equatorial band and transition zone (ztilde)
241                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
242!                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
243             
244                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
245                     ! values inside the equatorial band (ztilde as zstar)
246                     frq_rst_e3t(ji,jj) =  0.0_wp
247!                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
248                     tildemask(ji,jj) = 0._wp
249                  ELSE
250                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
251                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
252                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
253                        &                                          * 180._wp / 3.5_wp ) )
254!                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
255!                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
256!                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
257!                        &                                          * 180._wp / 3.5_wp ) )
258                     tildemask(ji,jj) = 0.5_wp * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
259                        &                                                 * 180._wp / 3.5_wp ) )
260                  ENDIF
261               END DO
262            END DO
263         ENDIF
264         !
265         IF ( ln_vvl_zstar_on_shelf ) THEN
266            zhmin = 50._wp
267            zhmax = 100._wp
268            DO jj = 1, jpj
269               DO ji = 1, jpi
270                  zwgt = 1._wp
271                  IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
272                     zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
273                  ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
274                     zwgt = 0._wp
275                  ENDIF
276                  frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt)
277                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt)
278               END DO
279            END DO
280         ENDIF
281         !
282         ztmp = MAXVAL( frq_rst_hdv(:,:) )
283         IF( lk_mpp )   CALL mpp_max( ztmp )                 ! max over the global domain
284         !
285         IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' )
286         !
287      ENDIF
288
289      IF( ln_vvl_layer ) THEN
290         IF ( ln_vvl_zstar_on_shelf ) THEN
291            zhmin = 100._wp
292            zhmax = 150._wp
293            DO jj = 1, jpj
294               DO ji = 1, jpi
295                  zwgt = 1._wp
296                  IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
297                     zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
298                  ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
299                     zwgt = 0._wp
300                  ENDIF
301                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt)
302               END DO
303            END DO
304         ENDIF
305         IF ( ln_vvl_zstar_at_eqtor ) THEN
306            DO jj = 1, jpj
307               DO ji = 1, jpi
308!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
309                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
310                     ! values outside the equatorial band and transition zone (ztilde)
311             
312                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
313                     ! values inside the equatorial band (ztilde as zstar)
314                     tildemask(ji,jj) = 0._wp
315                  ELSE
316                     tildemask(ji,jj) = 0.5_wp * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
317                        &                                                 * 180._wp / 3.5_wp ) )
318                  ENDIF
319               END DO
320            END DO
321         ENDIF
322      ENDIF
323      !
324      IF(lwxios) THEN
325! define variables in restart file when writing with XIOS
326         CALL iom_set_rstw_var_active('e3t_b')
327         CALL iom_set_rstw_var_active('e3t_n')
328         !                                           ! ----------------------- !
329         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
330            !                                        ! ----------------------- !
331            CALL iom_set_rstw_var_active('tilde_e3t_b')
332            CALL iom_set_rstw_var_active('tilde_e3t_n')
333         END IF
334         !                                           ! -------------!   
335         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
336            !                                        ! ------------ !
337            CALL iom_set_rstw_var_active('un_lf')
338            CALL iom_set_rstw_var_active('vn_lf')
339            CALL iom_set_rstw_var_active('hdivn_lf')
340         ENDIF
341         !
342      ENDIF
343      !
344   END SUBROUTINE dom_vvl_init
345
346
347   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
348      !!----------------------------------------------------------------------
349      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
350      !!                   
351      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
352      !!                 tranxt and dynspg routines
353      !!
354      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
355      !!               - z_tilde_case: after scale factor increment =
356      !!                                    high frequency part of horizontal divergence
357      !!                                  + retsoring towards the background grid
358      !!                                  + thickness difusion
359      !!                               Then repartition of ssh INCREMENT proportionnaly
360      !!                               to the "baroclinic" level thickness.
361      !!
362      !! ** Action  :  - hdivn_lf   : restoring towards full baroclinic divergence in z_tilde case
363      !!               - tilde_e3t_a: after increment of vertical scale factor
364      !!                              in z_tilde case
365      !!               - e3(t/u/v)_a
366      !!
367      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
368      !!----------------------------------------------------------------------
369      INTEGER, INTENT( in )           ::   kt      ! time step
370      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
371      !
372      LOGICAL                ::   ll_do_bclinic         ! local logical
373      INTEGER                ::   ji, jj, jk, jo        ! dummy loop indices
374      INTEGER                ::   ib, ib_bdy, ip, jp    !   "     "     "
375      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
376      INTEGER                ::   ncall
377      REAL(wp)               ::   z2dt  , z_tmin, z_tmax! local scalars               
378      REAL(wp)               ::   zalpha, zwgt          ! temporary scalars
379      REAL(wp)               ::   zdu, zdv, zramp, zmet
380      REAL(wp)               ::   ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv
381      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
382      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t, ztu, ztv
383      !!----------------------------------------------------------------------
384      !
385      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
386      !
387      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
388      !
389      IF( kt == nit000 ) THEN
390         IF(lwp) WRITE(numout,*)
391         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
392         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
393      ENDIF
394
395      zmet = 1._wp
396
397      ll_do_bclinic = .TRUE.
398      ncall = 1
399
400      IF( PRESENT(kcall) ) THEN
401         ! comment line below if tilda coordinate has to be computed at each call
402         IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE.
403         ncall = kcall 
404      ENDIF
405
406      IF( neuler == 0 .AND. kt == nit000 ) THEN
407         z2dt =  rdt
408      ELSE
409         z2dt = 2.0_wp * rdt
410      ENDIF
411
412      ! ******************************* !
413      ! After scale factors at t-points !
414      ! ******************************* !
415      !                                                ! --------------------------------------------- !
416      !                                                ! z_star coordinate and barotropic z-tilde part !
417      !                                                ! --------------------------------------------- !
418      !
419      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
420      DO jk = 1, jpkm1
421         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
422         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
423      END DO
424      !
425      IF ((ln_vvl_ztilde.OR.ln_vvl_layer).AND.(zmet==1._wp)) THEN
426         DO jk = 1, jpkm1
427            e3t_a(:,:,jk) = e3t_a(:,:,jk) - tilde_e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
428         END DO 
429      ENDIF
430
431      IF( (ln_vvl_ztilde.OR.ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
432         !                                                             ! ------baroclinic part------ !
433         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
434         un_td(:,:,:) = 0.0_wp        ! Transport corrections
435         vn_td(:,:,:) = 0.0_wp
436
437         zhdiv(:,:) = 0.
438         DO jk = 1, jpkm1
439            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
440         END DO
441         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) ) 
442
443         ze3t(:,:,:) = 0._wp
444         IF( ln_rnf ) THEN
445            CALL sbc_rnf_div( ze3t )          ! runoffs
446            DO jk=1,jpkm1
447               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ze3t(:,:,jk) * e3t_n(:,:,jk)
448            END DO   
449         ENDIF 
450
451         ! Thickness advection:
452         ! --------------------
453         ! Set advection velocities and source term
454         IF ( ln_vvl_ztilde ) THEN
455            IF ( ncall==1 ) THEN
456               zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
457               DO jk = 1, jpkm1
458                  ztu(:,:,jk) = un(:,:,jk) * e3u_n(:,:,jk) * e2u(:,:)
459                  ztv(:,:,jk) = vn(:,:,jk) * e3v_n(:,:,jk) * e1v(:,:)
460                  ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - (e3t_n(:,:,jk)-zmet*tilde_e3t_n(:,:,jk)) * zhdiv(:,:)
461               END DO
462               !
463                  un_lf(:,:,:,nn_filt_order)  =    un_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ztu(:,:,:)
464                  vn_lf(:,:,:,nn_filt_order)  =    vn_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ztv(:,:,:)
465               hdivn_lf(:,:,:,nn_filt_order)  = hdivn_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
466
467               DO jo = nn_filt_order-1,1,-1
468                     un_lf(:,:,:,jo)  =    un_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha *    un_lf(:,:,:,jo+1)
469                     vn_lf(:,:,:,jo)  =    vn_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha *    vn_lf(:,:,:,jo+1)
470                  hdivn_lf(:,:,:,jo)  = hdivn_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha * hdivn_lf(:,:,:,jo+1)
471               END DO
472            ENDIF
473            !
474            DO jk = 1, jpkm1
475               ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk,1)/e3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk)
476               ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk,1)/e3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk)
477               tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk,1) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
478            END DO
479            !
480         ELSEIF ( ln_vvl_layer ) THEN
481            !
482            DO jk = 1, jpkm1
483               ztu(:,:,jk) = un(:,:,jk)
484               ztv(:,:,jk) = vn(:,:,jk)
485            END DO
486            !
487         ENDIF
488         !
489         ! Block fluxes through small layers:
490!         DO jk=1,jpkm1
491!            DO ji = 1, jpi
492!               DO jj= 1, jpj
493!                  zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_n(ji,jj,jk) - hsmall) )
494!                  un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * e3u_n(ji,jj,jk) * e2u(ji,jj)
495!                  ztu(ji,jj,jk) = zmsku * ztu(ji,jj,jk)
496!                  IF ( ln_vvl_ztilde ) un_lf(ji,jj,jk) = zmsku * un_lf(ji,jj,jk)
497!                  !
498!                  zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_n(ji,jj,jk) - hsmall) )
499!                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * e3v_n(ji,jj,jk) * e1v(ji,jj)
500!                  ztv(ji,jj,jk) = zmskv * ztv(ji,jj,jk)
501!                  IF ( ln_vvl_ztilde ) vn_lf(ji,jj,jk) = zmskv * vn_lf(ji,jj,jk)
502!               END DO
503!            END DO
504!         END DO     
505         !
506         ! Do advection
507         IF     (ln_vvl_adv_fct) THEN
508            CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv )
509            !
510         ELSEIF (ln_vvl_adv_cn2) THEN
511            DO jk = 1, jpkm1
512               DO jj = 2, jpjm1
513                  DO ji = fs_2, fs_jpim1   ! vector opt.
514                     tilde_e3t_a(ji,jj,jk) =   tilde_e3t_a(ji,jj,jk) &
515                     & -(  e2u(ji,jj)*e3u_n(ji,jj,jk) * ztu(ji,jj,jk) - e2u(ji-1,jj  )*e3u_n(ji-1,jj  ,jk) * ztu(ji-1,jj  ,jk)       &
516                     &   + e1v(ji,jj)*e3v_n(ji,jj,jk) * ztv(ji,jj,jk) - e1v(ji  ,jj-1)*e3v_n(ji  ,jj-1,jk) * ztv(ji  ,jj-1,jk)  )    &
517                     & / ( e1t(ji,jj) * e2t(ji,jj) )
518                  END DO
519               END DO
520            END DO
521         ENDIF
522         !
523         ! Thickness anomaly diffusion:
524         ! ----------------------------
525         ztu(:,:,:) = 0.0_wp
526         ztv(:,:,:) = 0.0_wp
527
528         ze3t(:,:,1) = 0.e0
529         DO jj = 1, jpj
530            DO ji = 1, jpi
531               DO jk = 2, jpk
532                  ze3t(ji,jj,jk) =  ze3t(ji,jj,jk-1) + tilde_e3t_b(ji,jj,jk-1) * tmask(ji,jj,jk-1) 
533               END DO
534            END DO
535         END DO
536
537         IF ( ln_vvl_blp ) THEN  ! Bilaplacian
538            DO jk = 1, jpkm1
539               DO jj = 1, jpjm1                 ! First derivative (gradient)
540                  DO ji = 1, fs_jpim1   ! vector opt.                 
541                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
542                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) )
543                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
544                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) )
545                  END DO
546               END DO
547
548               DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
549                  DO ji = fs_2, fs_jpim1 ! vector opt.
550                     zht(ji,jj) =  rn_ahe3_blp * r1_e1e2t(ji,jj) * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
551                        &                                            + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
552                  END DO
553               END DO
554
555               ! Open boundary conditions:
556               IF ( ln_bdy ) THEN
557                  DO ib_bdy=1, nb_bdy
558                     DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
559                        ji = idx_bdy(ib_bdy)%nbi(ib,1)
560                        jj = idx_bdy(ib_bdy)%nbj(ib,1)
561                        zht(ji,jj) = 0._wp
562                     END DO
563                  END DO
564               END IF
565
566               CALL lbc_lnk( zht, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
567
568               DO jj = 1, jpjm1                 ! third derivative (gradient)
569                  DO ji = 1, fs_jpim1   ! vector opt.
570                     ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) )
571                     ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) * ( zht(ji  ,jj+1) - zht(ji,jj) )
572                  END DO
573               END DO
574            END DO
575         ENDIF
576
577         IF ( ln_vvl_lap ) THEN    ! Laplacian
578            DO jk = 1, jpkm1                    ! First derivative (gradient)
579               DO jj = 1, jpjm1
580                  DO ji = 1, fs_jpim1   ! vector opt.                 
581                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) &
582                         &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) )
583                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
584                         &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) )
585                     ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu
586                     ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv
587                  END DO
588               END DO
589            END DO
590         ENDIF 
591
592         ! divergence of diffusive fluxes
593         DO jk = 1, jpkm1
594            DO jj = 2, jpjm1
595               DO ji = fs_2, fs_jpim1   ! vector opt.
596                  un_td(ji,jj,jk) = un_td(ji,jj,jk) + ztu(ji,jj,jk+1) - ztu(ji,jj,jk  ) 
597                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) + ztv(ji,jj,jk+1) - ztv(ji,jj,jk  ) 
598                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk+1) - ztu(ji,jj,jk+1)    &
599                     &                                               +ztv(ji  ,jj-1,jk+1) - ztv(ji,jj,jk+1)    &
600                     &                                               -ztu(ji-1,jj  ,jk  ) + ztu(ji,jj,jk  )    &
601                     &                                               -ztv(ji  ,jj-1,jk  ) + ztv(ji,jj,jk  )    &
602                     &                                            ) * r1_e1e2t(ji,jj)
603               END DO
604            END DO
605         END DO
606
607         CALL lbc_lnk_multi( un_td, 'U', -1., vn_td, 'V', -1. )     !* local domain boundaries
608         !
609         CALL dom_vvl_ups_cor( kt, tilde_e3t_a, un_td, vn_td )
610
611!         IF ( ln_vvl_ztilde ) THEN
612!            ztu(:,:,:) = -un_lf(:,:,:)
613!            ztv(:,:,:) = -vn_lf(:,:,:)
614!            CALL dom_vvl_ups_cor( kt, tilde_e3t_a, ztu, ztv )
615!         ENDIF
616         !
617         ! Remove "external thickness" tendency:
618         DO jk = 1, jpkm1
619            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   (e3t_n(:,:,jk)-zmet*tilde_e3t_n(:,:,jk)) * zhdiv(:,:) * tmask(:,:,jk)
620         END DO 
621                   
622         ! Leapfrog time stepping
623         ! ~~~~~~~~~~~~~~~~~~~~~~
624         zramp = 1._wp
625         IF ((.NOT.ln_rstart).AND.ln_vvl_ramp) zramp = MIN(MAX( ((kt-nit000)*rdt)/(rn_day_ramp*rday),0._wp),1._wp)
626
627         DO jk=1,jpkm1
628            tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) &
629                               & * tildemask(:,:) * zramp
630         END DO
631
632         ! Ensure layer separation:
633         ! ~~~~~~~~~~~~~~~~~~~~~~~~
634         IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt )
635 
636         ! Boundary conditions:
637         ! ~~~~~~~~~~~~~~~~~~~~
638         IF ( ln_bdy ) THEN
639            DO ib_bdy = 1, nb_bdy
640               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
641!!               DO ib = 1, idx_bdy(ib_bdy)%nblen(1)
642                  ji = idx_bdy(ib_bdy)%nbi(ib,1)
643                  jj = idx_bdy(ib_bdy)%nbj(ib,1)
644                  zwgt = idx_bdy(ib_bdy)%nbw(ib,1)
645                  ip = bdytmask(ji+1,jj  ) - bdytmask(ji-1,jj  )
646                  jp = bdytmask(ji  ,jj+1) - bdytmask(ji  ,jj-1)
647                  DO jk = 1, jpkm1 
648                     tilde_e3t_a(ji,jj,jk) = 0.e0
649!!                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt)
650!!                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk)
651                  END DO
652               END DO
653            END DO
654         ENDIF
655
656         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )       
657       
658      ENDIF
659
660      IF( ln_vvl_ztilde.AND.( ncall==1))  THEN
661         zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
662         !
663         ! divergence of diffusive fluxes
664         DO jk = 1, jpkm1
665            DO jj = 2, jpjm1
666               DO ji = fs_2, fs_jpim1   ! vector opt.
667                  ze3t(ji,jj,jk) = (  un_td(ji,jj,jk) - un_td(ji-1,jj  ,jk) &
668                                 &  + vn_td(ji,jj,jk) - vn_td(ji  ,jj-1,jk) &
669                                 & ) / ( e1t(ji,jj) * e2t(ji,jj) )
670               END DO
671            END DO
672         END DO
673         CALL lbc_lnk( ze3t(:,:,:), 'T', 1. )
674
675         DO jo = nn_filt_order,1,-1
676            hdivn_lf(:,:,:,jo) =  hdivn_lf(:,:,:,jo)  + zalpha**(nn_filt_order-jo+1) * ze3t(:,:,:)
677         END DO
678      ENDIF
679
680      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
681      !                                             ! ---baroclinic part--------- !
682
683         IF ( (ncall==2).AND.(.NOT.ll_do_bclinic) ) THEN
684
685            DO jk = 1, jpkm1
686               ztu(:,:,jk) = (un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * e3u_n(:,:,jk) * e2u(:,:) * umask(:,:,jk)
687               ztv(:,:,jk) = (vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * e3v_n(:,:,jk) * e1v(:,:) * vmask(:,:,jk)
688               DO jj = 2, jpjm1
689                  DO ji = fs_2, fs_jpim1   ! vector opt.
690                     ze3t(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk) &
691                                    &  + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk) &
692                                    & ) / ( e1t(ji,jj) * e2t(ji,jj) )
693                  END DO
694               END DO
695            END DO
696            !
697            zhdiv(:,:) = 0.
698            DO jk = 1, jpkm1
699               zhdiv(:,:) = zhdiv(:,:) + ze3t(:,:,jk) * tmask(:,:,jk)
700            END DO
701            zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) ) 
702            !
703            DO jk = 1, jpkm1
704               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * (ze3t(:,:,jk) & 
705                                   & - zhdiv(:,:)*(e3t_n(:,:,jk)-zmet*tilde_e3t_n(:,:,jk))*tmask(:,:,jk))
706            END DO
707            CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
708         ENDIF
709         !
710         DO jk = 1, jpkm1
711            e3t_a(:,:,jk) = e3t_a(:,:,jk) + tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
712         END DO
713!         IF( ln_vvl_ztilde ) THEN !  Relax barotropic component:
714!            DO jk = 1, jpkm1
715!               e3t_a(:,:,jk) = e3t_a(:,:,jk)  &
716!                             & - z2dt * frq_rst_e3t(:,:) * (e3t_b(:,:,jk) - tilde_e3t_b(:,:,jk)   &
717!                             & - e3t_0(:,:,jk) * (ht_0(:,:) + sshb(:,:))/ (ht_0(:,:)*tmask(:,:,1) + 1._wp - tmask(:,:,1)))                 
718!            END DO
719!         ENDIF
720      ENDIF
721
722      IF( ln_vvl_dbg.AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN   ! - ML - test: control prints for debuging
723         !
724         zht(:,:) = 0.0_wp
725         DO jk = 1, jpkm1
726            zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
727         END DO
728         IF( lwp ) WRITE(numout, *) 'kt = ', kt
729         IF( lwp ) WRITE(numout, *) 'ncall = ', ncall
730         IF( lwp ) WRITE(numout, *) 'll_do_bclinic', ll_do_bclinic 
731         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
732            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 )
733            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
734            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
735         END IF
736         !
737         z_tmin = MINVAL( e3t_n(:,:,:)  )
738         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
739         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3t_n) =', z_tmin
740         IF ( z_tmin .LE. 0._wp ) THEN
741            IF( lk_mpp ) THEN
742               CALL mpp_minloc(e3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
743            ELSE
744               ijk_min = MINLOC( e3t_n(:,:,:)  )
745               ijk_min(1) = ijk_min(1) + nimpp - 1
746               ijk_min(2) = ijk_min(2) + njmpp - 1
747            ENDIF
748            IF (lwp) THEN
749               ji = ijk_min(1) ; jj = ijk_min(2) ; jk = ijk_min(3)
750               WRITE(numout, *) 'Negative scale factor, e3t_n =', z_tmin
751               WRITE(numout, *) 'at i, j, k=', ijk_min           
752               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
753            ENDIF
754         ENDIF         
755         !
756         z_tmin = MINVAL( e3u_n(:,:,:))
757         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
758         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3u_n) =', z_tmin
759         IF ( z_tmin .LE. 0._wp ) THEN
760            IF( lk_mpp ) THEN
761               CALL mpp_minloc(e3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
762            ELSE
763               ijk_min = MINLOC( e3u_n(:,:,:)  )
764               ijk_min(1) = ijk_min(1) + nimpp - 1
765               ijk_min(2) = ijk_min(2) + njmpp - 1
766            ENDIF
767            IF (lwp) THEN
768               WRITE(numout, *) 'Negative scale factor, e3u_n =', z_tmin
769               WRITE(numout, *) 'at i, j, k=', ijk_min           
770               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
771            ENDIF
772         ENDIF 
773         !
774         z_tmin = MINVAL( e3v_n(:,:,:) )
775         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
776         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3v_n) =', z_tmin
777         IF ( z_tmin .LE. 0._wp ) THEN
778            IF( lk_mpp ) THEN
779               CALL mpp_minloc(e3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
780            ELSE
781               ijk_min = MINLOC( e3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0   )
782               ijk_min(1) = ijk_min(1) + nimpp - 1
783               ijk_min(2) = ijk_min(2) + njmpp - 1
784            ENDIF
785            IF (lwp) THEN
786               WRITE(numout, *) 'Negative scale factor, e3v_n =', z_tmin
787               WRITE(numout, *) 'at i, j, k=', ijk_min           
788               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
789            ENDIF
790         ENDIF 
791         !
792         z_tmin = MINVAL( e3f_n(:,:,:))
793         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
794         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3f_n) =', z_tmin
795         IF ( z_tmin .LE. 0._wp ) THEN
796            IF( lk_mpp ) THEN
797               CALL mpp_minloc(e3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
798            ELSE
799               ijk_min = MINLOC( e3f_n(:,:,:) )
800               ijk_min(1) = ijk_min(1) + nimpp - 1
801               ijk_min(2) = ijk_min(2) + njmpp - 1
802            ENDIF
803            IF (lwp) THEN
804               WRITE(numout, *) 'Negative scale factor, e3f_n =', z_tmin
805               WRITE(numout, *) 'at i, j, k=', ijk_min           
806               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
807            ENDIF
808         ENDIF
809         !
810         zht(:,:) = 0.0_wp
811         DO jk = 1, jpkm1
812            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
813         END DO
814         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
815         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
816         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn  -SUM(e3t_n))) =', z_tmax
817         !
818         zht(:,:) = 0.0_wp
819         DO jk = 1, jpkm1
820            zht(:,:) = zht(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
821         END DO
822         zwu(:,:) = 0._wp
823         DO jj = 1, jpjm1
824            DO ji = 1, fs_jpim1   ! vector opt.
825               zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e1e2u(ji,jj)                             &
826                        &          * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj) )
827            END DO
828         END DO
829         CALL lbc_lnk( zwu(:,:), 'U', 1._wp )
830         z_tmax = MAXVAL( ssumask(:,:) * ssumask(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) )
831         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
832         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(e3u_n))) =', z_tmax
833         !
834         zht(:,:) = 0.0_wp
835         DO jk = 1, jpkm1
836            zht(:,:) = zht(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
837         END DO
838         zwv(:,:) = 0._wp
839         DO jj = 1, jpjm1
840            DO ji = 1, fs_jpim1   ! vector opt.
841               zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e1e2v(ji,jj)                             &
842                        &          * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1) )
843            END DO
844         END DO
845         CALL lbc_lnk( zwv(:,:), 'V', 1._wp )
846         z_tmax = MAXVAL( ssvmask(:,:) * ssvmask(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) )
847         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
848         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(e3v_n))) =', z_tmax
849         !
850         zht(:,:) = 0.0_wp
851         DO jk = 1, jpkm1
852            DO jj = 1, jpjm1
853               zht(:,jj) = zht(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
854            END DO
855         END DO
856         zwu(:,:) = 0._wp
857         DO jj = 1, jpjm1
858            DO ji = 1, fs_jpim1   ! vector opt.
859               zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e1e2f(ji,jj)  &
860                        &          * (  e1e2t(ji  ,jj) * sshn(ji  ,jj) + e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1) &
861                        &             + e1e2t(ji+1,jj) * sshn(ji+1,jj) + e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) )
862            END DO
863         END DO
864         CALL lbc_lnk( zht(:,:), 'F', 1._wp )
865         CALL lbc_lnk( zwu(:,:), 'F', 1._wp )
866         z_tmax = MAXVAL( fmask(:,:,1) * fmask(:,:,1) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) )
867         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
868         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(e3f_n))) =', z_tmax
869         !
870      END IF
871
872      ! *********************************** !
873      ! After scale factors at u- v- points !
874      ! *********************************** !
875
876      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
877      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
878
879      ! *********************************** !
880      ! After depths at u- v points         !
881      ! *********************************** !
882
883      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
884      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
885      DO jk = 2, jpkm1
886         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
887         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
888      END DO
889      !                                        ! Inverse of the local depth
890!!gm BUG ?  don't understand the use of umask_i here .....
891      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
892      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
893      !
894      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
895      !
896   END SUBROUTINE dom_vvl_sf_nxt
897
898
899   SUBROUTINE dom_vvl_sf_swp( kt )
900      !!----------------------------------------------------------------------
901      !!                ***  ROUTINE dom_vvl_sf_swp  ***
902      !!                   
903      !! ** Purpose :  compute time filter and swap of scale factors
904      !!               compute all depths and related variables for next time step
905      !!               write outputs and restart file
906      !!
907      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
908      !!               - reconstruct scale factor at other grid points (interpolate)
909      !!               - recompute depths and water height fields
910      !!
911      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
912      !!               - Recompute:
913      !!                    e3(u/v)_b       
914      !!                    e3w_n           
915      !!                    e3(u/v)w_b     
916      !!                    e3(u/v)w_n     
917      !!                    gdept_n, gdepw_n  and gde3w_n
918      !!                    h(u/v) and h(u/v)r
919      !!
920      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
921      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
922      !!----------------------------------------------------------------------
923      INTEGER, INTENT( in ) ::   kt   ! time step
924      !
925      INTEGER  ::   ji, jj, jk   ! dummy loop indices
926      REAL(wp) ::   zcoef        ! local scalar
927      !!----------------------------------------------------------------------
928      !
929      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
930      !
931      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
932      !
933      IF( kt == nit000 )   THEN
934         IF(lwp) WRITE(numout,*)
935         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
936         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
937      ENDIF
938      !
939      ! Time filter and swap of scale factors
940      ! =====================================
941      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
942      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
943         IF( neuler == 0 .AND. kt == nit000 ) THEN
944            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
945         ELSE
946            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
947            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
948         ENDIF
949         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
950      ENDIF
951      gdept_b(:,:,:) = gdept_n(:,:,:)
952      gdepw_b(:,:,:) = gdepw_n(:,:,:)
953
954      e3t_n(:,:,:) = e3t_a(:,:,:)
955      e3u_n(:,:,:) = e3u_a(:,:,:)
956      e3v_n(:,:,:) = e3v_a(:,:,:)
957
958      ! Compute all missing vertical scale factor and depths
959      ! ====================================================
960      ! Horizontal scale factor interpolations
961      ! --------------------------------------
962      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
963      ! - JC - hu_b, hv_b, hur_b, hvr_b also
964     
965      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F'  )
966     
967      ! Vertical scale factor interpolations
968      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
969      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
970      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
971      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
972      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
973      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
974
975      ! t- and w- points depth (set the isf depth as it is in the initial step)
976      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
977      gdepw_n(:,:,1) = 0.0_wp
978      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
979      DO jk = 2, jpk
980         DO jj = 1,jpj
981            DO ji = 1,jpi
982              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
983                                                                 ! 1 for jk = mikt
984               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
985               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
986               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
987                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
988               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
989            END DO
990         END DO
991      END DO
992
993      ! Local depth and Inverse of the local depth of the water
994      ! -------------------------------------------------------
995      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
996      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
997      !
998      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
999      DO jk = 2, jpkm1
1000         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
1001      END DO
1002
1003      ! write restart file
1004      ! ==================
1005      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
1006      !
1007      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
1008      !
1009   END SUBROUTINE dom_vvl_sf_swp
1010
1011
1012   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
1013      !!---------------------------------------------------------------------
1014      !!                  ***  ROUTINE dom_vvl__interpol  ***
1015      !!
1016      !! ** Purpose :   interpolate scale factors from one grid point to another
1017      !!
1018      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
1019      !!                - horizontal interpolation: grid cell surface averaging
1020      !!                - vertical interpolation: simple averaging
1021      !!----------------------------------------------------------------------
1022      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
1023      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
1024      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
1025      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1026      !
1027      INTEGER ::   ji, jj, jk, jkbot                                ! dummy loop indices
1028      INTEGER ::   nmet                                             ! horizontal interpolation method
1029      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
1030      REAL(wp) ::  ztap, zsmall                                     ! Parameters defining minimum thicknesses UVF-points
1031      REAL(wp) ::  zmin
1032      REAL(wp) ::  zdo, zup                                         ! Lower and upper interfaces depths anomalies
1033      REAL(wp), DIMENSION(jpi,jpj) :: zs                            ! Surface interface depth anomaly
1034      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw                        ! Interface depth anomaly
1035      !!----------------------------------------------------------------------
1036      !
1037!     nmet = 0  ! Original method (Surely wrong)
1038     nmet = 2  ! Interface interpolation
1039!      nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment
1040!     Note that we kept surface weighted interpolation for barotropic increment to be compliant
1041!     with what is done in surface pressure module.
1042      !
1043      IF(ln_wd_il) THEN
1044        zlnwd = 1.0_wp
1045      ELSE
1046        zlnwd = 0.0_wp
1047      END IF
1048      !
1049      ztap   = 0._wp   ! Minimum fraction of T-point thickness at cell interfaces
1050      zsmall = 1.e-8_wp ! Minimum thickness at U or V points (m)
1051      !
1052      IF ( (nmet==1).OR.(nmet==2) ) THEN
1053         SELECT CASE ( pout )
1054            !
1055         CASE( 'U', 'V', 'F' )
1056            ! Compute interface depth anomaly at T-points
1057            !
1058            zw(:,:,:) =  0._wp   
1059            !
1060            DO jk=2,jpk
1061               zw(:,:,jk) = zw(:,:,jk-1) + pe3_in(:,:,jk-1)*tmask(:,:,jk-1)
1062            END DO 
1063            ! Interface depth anomalies:
1064            DO jk=1,jpkm1
1065               zw(:,:,jk) = zw(:,:,jk) - zw(:,:,jpk) + ht_0(:,:)
1066            END DO
1067            zw(:,:,jpk) = ht_0(:,:)
1068            !
1069            IF (nmet==2) THEN        ! Consider "internal" interfaces only
1070               zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh)
1071               !
1072               DO jj = 1, jpj
1073                  DO ji = 1, jpi
1074                     DO jk=1,jpk
1075                        zw(ji,jj,jk) = (zw(ji,jj,jk) + zs(ji,jj))                                          &
1076                                     & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1))  &
1077                                     & * tmask(ji,jj,jk)
1078                     END DO
1079                  END DO
1080               END DO
1081            ENDIF
1082            zw(:,:,:) = (zw(:,:,:) - gdepw_0(:,:,:))*tmask(:,:,:)
1083            !
1084         END SELECT
1085      END IF
1086      !
1087      pe3_out(:,:,:) = 0.0_wp
1088      !
1089      SELECT CASE ( pout )    !==  type of interpolation  ==!
1090         !
1091      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
1092         IF (nmet==0) THEN
1093            DO jk = 1, jpk
1094               DO jj = 1, jpjm1
1095                  DO ji = 1, fs_jpim1   ! vector opt.
1096                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
1097                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
1098                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
1099                  END DO
1100               END DO
1101            END DO
1102         ELSE
1103            DO jj = 1, jpjm1
1104               DO ji = 1, fs_jpim1   ! vector opt.
1105                  ! Correction at last level:
1106                  jkbot = mbku(ji,jj)
1107                  zdo = 0._wp
1108                  DO jk=jkbot,1,-1
1109                     zup = 0.5_wp * ( e1e2t(ji  ,jj)*zw(ji  ,jj,jk) &
1110                         &          + e1e2t(ji+1,jj)*zw(ji+1,jj,jk) ) * r1_e1e2u(ji,jj)
1111                     !
1112                     ! If there is a step, taper bottom interface:
1113!                     IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(jk==jkbot)) THEN
1114!                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
1115!                           zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))
1116!                        ELSE
1117!                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
1118!                        ENDIF
1119!                        zup = MIN(zup, zdo-zmin)
1120!                     ENDIF
1121                     zup = MIN(zup, zdo+e3u_0(ji,jj,jk)-zsmall)
1122                     pe3_out(ji,jj,jk) = (zdo - zup) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
1123                     zdo = zup
1124                  END DO
1125               END DO
1126            END DO
1127         END IF
1128         !
1129         IF (nmet==2) THEN           ! Spread sea level anomaly
1130            DO jj = 1, jpjm1
1131               DO ji = 1, fs_jpim1   ! vector opt.
1132                  DO jk=1,jpk
1133                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                   &
1134                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )               & 
1135                           &               / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )              &
1136                           &               * 0.5_wp * r1_e1e2u(ji,jj)                              &
1137                           &               * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )        &
1138                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji+1,jj)*zs(ji+1,jj) )       
1139                  END DO
1140               END DO
1141            END DO
1142            !
1143         ENDIF
1144         !
1145         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
1146         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
1147         !
1148      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
1149         IF (nmet==0) THEN
1150            DO jk = 1, jpk
1151               DO jj = 1, jpjm1
1152                  DO ji = 1, fs_jpim1   ! vector opt.
1153                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
1154                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
1155                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
1156                  END DO
1157               END DO
1158            END DO
1159         ELSE
1160            DO jj = 1, jpjm1
1161               DO ji = 1, fs_jpim1   ! vector opt.
1162                  ! Correction at last level:
1163                  jkbot = mbkv(ji,jj)
1164                  zdo = 0._wp
1165                  DO jk=jkbot,1,-1
1166                     zup = 0.5_wp * ( e1e2t(ji,jj  ) * zw(ji,jj  ,jk) & 
1167                         &          + e1e2t(ji,jj+1) * zw(ji,jj+1,jk) ) * r1_e1e2v(ji,jj)
1168                     !
1169                     ! If there is a step, taper bottom interface:
1170!                     IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN
1171!                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
1172!                           zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))
1173!                        ELSE
1174!                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
1175!                        ENDIF
1176!                        zup = MIN(zup, zdo-zmin)
1177!                     ENDIF
1178                     zup = MIN(zup, zdo + e3v_0(ji,jj,jk) - zsmall)
1179                     pe3_out(ji,jj,jk) = (zdo - zup) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
1180                     zdo = zup
1181                  END DO
1182               END DO
1183            END DO
1184         END IF
1185         !
1186         IF (nmet==2) THEN           ! Spread sea level anomaly
1187            DO jj = 1, jpjm1
1188               DO ji = 1, fs_jpim1   ! vector opt.
1189                  DO jk=1,jpk
1190                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                                       &
1191                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )                                   & 
1192                           &               / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )                                  &
1193                           &               * 0.5_wp * r1_e1e2v(ji,jj)                                                  &
1194                                           * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )                            &
1195                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji,jj+1)*zs(ji,jj+1) )       
1196                  END DO
1197               END DO
1198            END DO
1199            !
1200         ENDIF
1201         !
1202         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
1203         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
1204         !
1205      CASE( 'F' )                   !* from T-point to F-point : hor. surface weighted mean
1206         IF (nmet==0) THEN
1207            DO jk=1,jpk
1208               DO jj = 1, jpjm1
1209                  DO ji = 1, fs_jpim1   ! vector opt.
1210                     pe3_out(ji,jj,jk) = 0.25_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
1211                           &                     *    r1_e1e2f(ji,jj)                                                  &
1212                           &                     * (  e1e2t(ji  ,jj  ) * ( pe3_in(ji  ,jj  ,jk)-e3t_0(ji  ,jj  ,jk) )  & 
1213                           &                        + e1e2t(ji  ,jj+1) * ( pe3_in(ji  ,jj+1,jk)-e3t_0(ji  ,jj+1,jk) )  &
1214                           &                        + e1e2t(ji+1,jj  ) * ( pe3_in(ji+1,jj  ,jk)-e3t_0(ji+1,jj  ,jk) )  & 
1215                           &                        + e1e2t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ) )                                                 
1216                  END DO
1217               END DO
1218            END DO
1219         ELSE
1220            DO jj = 1, jpjm1
1221               DO ji = 1, fs_jpim1   ! vector opt.
1222                  ! bottom correction:
1223                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 
1224                  zdo = 0._wp
1225                  DO jk=jkbot,1,-1
1226                     zup =  0.25_wp * (  e1e2t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  & 
1227                           &           + e1e2t(ji+1,jj  ) * zw(ji+1,jj  ,jk)  & 
1228                           &           + e1e2t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  & 
1229                           &           + e1e2t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  ) *    r1_e1e2f(ji,jj)
1230                     !
1231                     ! If there is a step, taper bottom interface:
1232!                     IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN
1233!                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
1234!                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN
1235!                              zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))
1236!                           ELSE
1237!                              zmin = ztap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk))
1238!                           ENDIF
1239!                        ELSE
1240!                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
1241!                              zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk))
1242!                           ELSE
1243!                              zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk))
1244!                           ENDIF
1245!                        ENDIF
1246!                        zup = MIN(zup, zdo-zmin)
1247!                     ENDIF
1248                     zup = MIN(zup, zdo+e3f_0(ji,jj,jk)-zsmall)
1249                     !
1250                     pe3_out(ji,jj,jk) = ( zdo - zup ) & 
1251                                      & *( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )
1252                     zdo = zup
1253                  END DO
1254               END DO
1255            END DO
1256         END IF
1257         !
1258         IF (nmet==2) THEN           ! Spread sea level anomaly
1259            !
1260            DO jj = 1, jpjm1
1261               DO ji = 1, fs_jpim1   ! vector opt.
1262                  DO jk=1,jpk
1263                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                           &
1264                           &               + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) )                       & 
1265                           &               / ( hf_0(ji,jj) + 1._wp - umask(ji,jj,1)*umask(ji,jj+1,1) )     &
1266                           &               * 0.25_wp * r1_e1e2f(ji,jj)                                     & 
1267                           &               * ( umask(ji,jj,jk)*umask(ji,jj+1,jk)*(1.0_wp - zlnwd) + zlnwd )&
1268                           &               * ( e1e2t(ji  ,jj)*zs(ji  ,jj) + e1e2t(ji  ,jj+1)*zs(ji  ,jj+1) &
1269                           &                  +e1e2t(ji+1,jj)*zs(ji+1,jj) + e1e2t(ji+1,jj+1)*zs(ji+1,jj+1) )               
1270                  END DO
1271               END DO
1272            END DO
1273         END IF
1274         !
1275         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
1276         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
1277         !
1278      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
1279         !
1280         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
1281         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
1282!!gm BUG? use here wmask in case of ISF ?  to be checked
1283         DO jk = 2, jpk
1284            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
1285               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
1286               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
1287               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
1288         END DO
1289         !
1290      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
1291         !
1292         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
1293         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1294!!gm BUG? use here wumask in case of ISF ?  to be checked
1295         DO jk = 2, jpk
1296            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
1297               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
1298               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
1299               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
1300         END DO
1301         !
1302      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
1303         !
1304         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
1305         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1306!!gm BUG? use here wvmask in case of ISF ?  to be checked
1307         DO jk = 2, jpk
1308            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
1309               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
1310               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
1311               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
1312         END DO
1313      END SELECT
1314      !
1315   END SUBROUTINE dom_vvl_interpol
1316
1317
1318   SUBROUTINE dom_vvl_rst( kt, cdrw )
1319      !!---------------------------------------------------------------------
1320      !!                   ***  ROUTINE dom_vvl_rst  ***
1321      !!                     
1322      !! ** Purpose :   Read or write VVL file in restart file
1323      !!
1324      !! ** Method  :   use of IOM library
1325      !!                if the restart does not contain vertical scale factors,
1326      !!                they are set to the _0 values
1327      !!                if the restart does not contain vertical scale factors increments (z_tilde),
1328      !!                they are set to 0.
1329      !!----------------------------------------------------------------------
1330      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1331      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1332      !
1333      INTEGER ::   ji, jj, jk
1334      INTEGER ::   id1, id2, id3, id4, id5, id6, id7     ! local integers
1335      !!----------------------------------------------------------------------
1336      !
1337      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1338         !                                   ! ===============
1339         IF( ln_rstart ) THEN                   !* Read the restart file
1340            CALL rst_read_open                  !  open the restart file if necessary
1341            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
1342            !
1343            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
1344            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
1345            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
1346            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
1347            id5 = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
1348            id6 = iom_varid( numror, 'un_lf', ldstop = .FALSE. )
1349            id7 = iom_varid( numror, 'vn_lf', ldstop = .FALSE. )
1350            !                             ! --------- !
1351            !                             ! all cases !
1352            !                             ! --------- !
1353            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
1354               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
1355               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
1356               ! needed to restart if land processor not computed
1357               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
1358               WHERE ( tmask(:,:,:) == 0.0_wp ) 
1359                  e3t_n(:,:,:) = e3t_0(:,:,:)
1360                  e3t_b(:,:,:) = e3t_0(:,:,:)
1361               END WHERE
1362               IF( neuler == 0 ) THEN
1363                  e3t_b(:,:,:) = e3t_n(:,:,:)
1364               ENDIF
1365            ELSE IF( id1 > 0 ) THEN
1366               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
1367               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
1368               IF(lwp) write(numout,*) 'neuler is forced to 0'
1369               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
1370               e3t_n(:,:,:) = e3t_b(:,:,:)
1371               neuler = 0
1372            ELSE IF( id2 > 0 ) THEN
1373               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
1374               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
1375               IF(lwp) write(numout,*) 'neuler is forced to 0'
1376               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
1377               e3t_b(:,:,:) = e3t_n(:,:,:)
1378               neuler = 0
1379            ELSE
1380               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
1381               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1382               IF(lwp) write(numout,*) 'neuler is forced to 0'
1383               DO jk = 1, jpk
1384                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
1385                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1386                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
1387               END DO
1388               e3t_b(:,:,:) = e3t_n(:,:,:)
1389               neuler = 0
1390            ENDIF
1391            !                             ! ----------- !
1392            IF( ln_vvl_zstar ) THEN       ! z_star case !
1393               !                          ! ----------- !
1394               IF( MIN( id3, id4 ) > 0 ) THEN
1395                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
1396               ENDIF
1397               !                          ! ----------------------- !
1398            ELSE                          ! z_tilde and layer cases !
1399               !                          ! ----------------------- !
1400               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
1401                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
1402                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
1403               ELSE                            ! one at least array is missing
1404                  tilde_e3t_b(:,:,:) = 0.0_wp
1405                  tilde_e3t_n(:,:,:) = 0.0_wp
1406               ENDIF
1407               !                          ! ------------ !
1408               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
1409                  !                       ! ------------ !
1410                  IF( MIN(id5, id6, id7) > 0 ) THEN  ! required arrays exist
1411                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:,1), ldxios = lrxios )
1412                     CALL iom_get( numror, jpdom_autoglo, 'un_lf', un_lf(:,:,:,1), ldxios = lrxios )
1413                     CALL iom_get( numror, jpdom_autoglo, 'vn_lf', vn_lf(:,:,:,1), ldxios = lrxios )
1414                  ELSE                ! array is missing
1415                     hdivn_lf(:,:,:,:) = 0.0_wp
1416                     un_lf(:,:,:,:) = 0.0_wp
1417                     vn_lf(:,:,:,:) = 0.0_wp
1418                  ENDIF
1419               ENDIF
1420            ENDIF
1421            !
1422         ELSE                                   !* Initialize at "rest"
1423            !
1424
1425            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
1426               !
1427               IF( cn_cfg == 'wad' ) THEN
1428                  ! Wetting and drying test case
1429                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
1430                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
1431                  sshn (:,:)     = sshb(:,:)
1432                  un   (:,:,:)   = ub  (:,:,:)
1433                  vn   (:,:,:)   = vb  (:,:,:)
1434               ELSE
1435                  ! if not test case
1436                  sshn(:,:) = -ssh_ref
1437                  sshb(:,:) = -ssh_ref
1438
1439                  DO jj = 1, jpj
1440                     DO ji = 1, jpi
1441                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
1442
1443                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1444                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1445                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1446                        ENDIF
1447                     ENDDO
1448                  ENDDO
1449               ENDIF !If test case else
1450
1451               ! Adjust vertical metrics for all wad
1452               DO jk = 1, jpk
1453                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
1454                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1455                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
1456               END DO
1457               e3t_b(:,:,:) = e3t_n(:,:,:)
1458
1459               DO ji = 1, jpi
1460                  DO jj = 1, jpj
1461                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
1462                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
1463                     ENDIF
1464                  END DO
1465               END DO 
1466               !
1467            ELSE
1468               !
1469               ! Just to read set ssh in fact, called latter once vertical grid
1470               ! is set up:
1471!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
1472!               !
1473!               DO jk=1,jpk
1474!                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
1475!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
1476!               END DO
1477!               e3t_n(:,:,:) = e3t_b(:,:,:)
1478                sshn(:,:)=0._wp
1479                e3t_n(:,:,:)=e3t_0(:,:,:)
1480                e3t_b(:,:,:)=e3t_0(:,:,:)
1481               !
1482            END IF           ! end of ll_wd edits
1483
1484            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1485               tilde_e3t_b(:,:,:) = 0._wp
1486               tilde_e3t_n(:,:,:) = 0._wp
1487               IF( ln_vvl_ztilde ) THEN
1488                  hdivn_lf(:,:,:,:) = 0._wp 
1489                  un_lf(:,:,:,:) = 0._wp 
1490                  vn_lf(:,:,:,:) = 0._wp 
1491               ENDIF
1492            END IF
1493         ENDIF
1494         !
1495      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1496         !                                   ! ===================
1497         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1498         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1499         !                                           ! --------- !
1500         !                                           ! all cases !
1501         !                                           ! --------- !
1502         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
1503         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
1504         !                                           ! ----------------------- !
1505         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
1506            !                                        ! ----------------------- !
1507            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
1508            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
1509         END IF
1510         !                                           ! -------------!   
1511         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
1512            !                                        ! ------------ !
1513            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:,1), ldxios = lwxios)
1514            CALL iom_rstput( kt, nitrst, numrow, 'un_lf', un_lf(:,:,:,1), ldxios = lwxios)
1515            CALL iom_rstput( kt, nitrst, numrow, 'vn_lf', vn_lf(:,:,:,1), ldxios = lwxios)
1516         ENDIF
1517         !
1518         IF( lwxios ) CALL iom_swap(      cxios_context          )
1519      ENDIF
1520      !
1521   END SUBROUTINE dom_vvl_rst
1522
1523
1524   SUBROUTINE dom_vvl_ctl
1525      !!---------------------------------------------------------------------
1526      !!                  ***  ROUTINE dom_vvl_ctl  ***
1527      !!               
1528      !! ** Purpose :   Control the consistency between namelist options
1529      !!                for vertical coordinate
1530      !!----------------------------------------------------------------------
1531      INTEGER ::   ioptio, ios
1532
1533      NAMELIST/nam_vvl/ ln_vvl_zstar               , ln_vvl_ztilde                       , &
1534                      & ln_vvl_layer               , ln_vvl_ztilde_as_zstar              , &
1535                      & ln_vvl_zstar_at_eqtor      , ln_vvl_zstar_on_shelf               , &
1536                      & ln_vvl_adv_cn2             , ln_vvl_adv_fct                      , &
1537                      & ln_vvl_lap                 , ln_vvl_blp                          , &
1538                      & rn_ahe3_lap                , rn_ahe3_blp                         , &
1539                      & rn_rst_e3t                 , rn_lf_cutoff                        , &
1540                      & ln_vvl_regrid                                                    , &
1541                      & ln_vvl_ramp                , rn_day_ramp                         , &
1542                      & ln_vvl_dbg   ! not yet implemented: ln_vvl_kepe
1543      !!---------------------------------------------------------------------- 
1544      !
1545      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
1546      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1547901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
1548      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
1549      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1550902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
1551      IF(lwm) WRITE ( numond, nam_vvl )
1552      !
1553      IF(lwp) THEN                    ! Namelist print
1554         WRITE(numout,*)
1555         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1556         WRITE(numout,*) '~~~~~~~~~~~'
1557         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
1558         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1559         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1560         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
1561         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1562         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
1563         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
1564         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1565         WRITE(numout,*) '      ztilde on shelves          ln_vvl_zstar_on_shelf  = ', ln_vvl_zstar_on_shelf
1566         WRITE(numout,*) '           Namelist nam_vvl : thickness advection scheme'
1567         WRITE(numout,*) '              2nd order                  ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2
1568         WRITE(numout,*) '              2nd order FCT              ln_vvl_adv_fct = ', ln_vvl_adv_fct                   
1569         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion scheme'
1570         WRITE(numout,*) '              Laplacian                  ln_vvl_lap     = ', ln_vvl_lap
1571         WRITE(numout,*) '              Bilaplacian                ln_vvl_blp     = ', ln_vvl_blp
1572         WRITE(numout,*) '              Laplacian   coefficient    rn_ahe3_lap    = ', rn_ahe3_lap
1573         WRITE(numout,*) '              Bilaplacian coefficient    rn_ahe3_blp    = ', rn_ahe3_blp
1574         WRITE(numout,*) '           Namelist nam_vvl : layers regriding'
1575         WRITE(numout,*) '              ln_vvl_regrid                             = ', ln_vvl_regrid
1576         WRITE(numout,*) '           Namelist nam_vvl : linear ramp at startup'
1577         WRITE(numout,*) '              ln_vvl_ramp                               = ', ln_vvl_ramp
1578         WRITE(numout,*) '              rn_day_ramp                               = ', rn_day_ramp
1579         IF( ln_vvl_ztilde_as_zstar ) THEN
1580            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
1581            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
1582            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
1583            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
1584            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1585            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
1586         ELSE
1587            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
1588            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
1589            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
1590            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
1591         ENDIF
1592         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
1593         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
1594      ENDIF
1595      !
1596      IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN
1597         ioptio = 0                      ! Choose one advection scheme at most
1598         IF( ln_vvl_adv_cn2         )        ioptio = ioptio + 1
1599         IF( ln_vvl_adv_fct         )        ioptio = ioptio + 1
1600         IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' )
1601      ENDIF
1602      !
1603      ioptio = 0                      ! Parameter control
1604      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1605      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1606      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1607      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1608      !
1609      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1610      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1611      !
1612      IF(lwp) THEN                   ! Print the choice
1613         WRITE(numout,*)
1614         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1615         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1616         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1617         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1618      ENDIF
1619      !
1620      ! Use of "shelf horizon depths" should be allowed with s-z coordinates, but we restrict it to zco and zps
1621      ! for the time being
1622      IF ( ln_sco ) THEN
1623        ll_shorizd=.FALSE.
1624      ELSE
1625        ll_shorizd=.TRUE.
1626      ENDIF
1627      !
1628#if defined key_agrif
1629      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1630#endif
1631      !
1632   END SUBROUTINE dom_vvl_ctl
1633
1634   SUBROUTINE dom_vvl_regrid( kt )
1635      !!----------------------------------------------------------------------
1636      !!                ***  ROUTINE dom_vvl_regrid  ***
1637      !!                   
1638      !! ** Purpose :  Ensure "well-behaved" vertical grid
1639      !!
1640      !! ** Method  :  More or less adapted from references below.
1641      !!regrid
1642      !! ** Action  :  Ensure that thickness are above a given value, spaced enough
1643      !!               and revert to Eulerian coordinates near the bottom.     
1644      !!
1645      !! References :  Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction
1646      !!               with a Model Combining Terrain-following and Isentropic
1647      !!               coordinates. Part I: Model Description. Monthly Weather Rev.,
1648      !!               121, 1770-1785.
1649      !!               Toy, M., 2011: Incorporating Condensational Heating into a
1650      !!               Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic-
1651      !!               Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954.
1652      !!----------------------------------------------------------------------
1653      !! * Arguments
1654      INTEGER, INTENT( in )               :: kt         ! time step
1655
1656      !! * Local declarations
1657      INTEGER                             :: ji, jj, jk ! dummy loop indices
1658      LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond
1659      LOGICAL                             :: ll_zdiff_cond, ll_blpdiff_cond
1660      INTEGER                             :: jkbot
1661      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd
1662      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj, dzmin_int, dzmin_surf
1663      REAL(wp)                            :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn
1664      REAL(wp)                            :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot
1665      REAL(wp)                            :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim
1666      REAL(wp), DIMENSION(jpi,jpj)        :: zdw, zwu, zwv
1667      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: zwdw, zwdw_b
1668      !!----------------------------------------------------------------------
1669
1670      IF( ln_timing )  CALL timing_start('dom_vvl_regrid')
1671      !
1672      !
1673      ! Some user defined parameters below:
1674      ll_chk_bot2top   = .TRUE.
1675      ll_chk_top2bot   = .TRUE.
1676      dzmin_int  = 1.0_wp   ! Absolute minimum depth in the interior (in meters)
1677      dzmin_surf = 1.0_wp   ! Absolute minimum depth at the surface (in meters)
1678      zfrch_stp  = 5._wp   ! Maximum fractionnal thickness change in one time step (<= 1.)
1679      zfrch_rel  = 0.4_wp   ! Maximum relative thickness change in the vertical (<= 1.)
1680      zfrac_bot  = 0.05_wp  ! Fraction of bottom level allowed to change
1681      zscal_bot  = 2.0_wp   ! Depth lengthscale
1682      ll_zdiff_cond    = .TRUE.  ! Conditionnal vertical diffusion of interfaces
1683         zvdiff  = 0.2_wp   ! m
1684         zvlim   = 0.5_wp   ! max d2h/dh
1685      ll_lapdiff_cond  = .TRUE.  ! Conditionnal Laplacian diffusion of interfaces
1686         zhdiff  = 0.01_wp  ! ad.
1687         zhlim   = 0.03_wp  ! ad. max lap(z)*e1
1688      ll_blpdiff_cond  = .TRUE.  ! Conditionnal Bilaplacian diffusion of interfaces
1689         zhdiff2 = 0.2_wp   ! ad.
1690         zhlim2  = 0.01_wp  ! ad. max bilap(z)*e1**3
1691      ! ---------------------------------------------------------------------------------------
1692      !
1693      ! Set arrays determining maximum vertical displacement at the bottom:
1694      !--------------------------------------------------------------------
1695      IF ( kt==nit000 ) THEN
1696         DO jj = 2, jpjm1
1697            DO ji = 2, jpim1
1698               jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1))
1699               jk = MIN(jk,mbkt(ji-1,jj-1), mbkt(ji-1,jj+1), mbkt(ji+1,jj+1), mbkt(ji+1,jj-1))
1700               i_int_bot(ji,jj) = jk
1701            END DO
1702         END DO
1703         dsm(:,:) = REAL( i_int_bot(:,:), wp )   ;   CALL lbc_lnk(dsm(:,:),'T',1.)
1704         i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 )
1705
1706         DO jj = 2, jpjm1
1707            DO ji = 2, jpim1
1708               zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji  ,jj,1), &
1709                          &     ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), &
1710                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj  ,1), &
1711                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1)  )
1712                zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall )
1713            END DO
1714         END DO
1715         CALL lbc_lnk( zdw(:,:), 'T', 1. )
1716
1717         DO jj = 2, jpjm1
1718            DO ji = 2, jpim1
1719               dsm(ji,jj) = 1._wp/16._wp * (            zdw(ji-1,jj-1) + zdw(ji+1,jj-1)      &
1720                  &                           +         zdw(ji-1,jj+1) + zdw(ji+1,jj+1)      &
1721                  &                           + 2._wp*( zdw(ji  ,jj-1) + zdw(ji-1,jj  )      &
1722                  &                           +         zdw(ji+1,jj  ) + zdw(ji  ,jj+1) )    &
1723                  &                           + 4._wp*  zdw(ji  ,jj  )                       )
1724            END DO
1725         END DO         
1726
1727         CALL lbc_lnk( dsm(:,:), 'T', 1. )   
1728     
1729         IF (ln_zps) THEN
1730            DO jj = 1, jpj
1731               DO ji = 1, jpi
1732                  jk = i_int_bot(ji,jj)
1733                  hsm(ji,jj) = zfrac_bot * e3w_1d(jk)
1734                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.05_wp*ht_0(ji,jj))
1735               END DO
1736            END DO
1737         ELSE
1738            DO jj = 1, jpj
1739               DO ji = 1, jpi
1740                  jk = i_int_bot(ji,jj)
1741                  hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk)
1742                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.05_wp*ht_0(ji,jj))
1743               END DO
1744            END DO
1745         ENDIF
1746      END IF
1747
1748      ! Provisionnal interface depths:
1749      !-------------------------------
1750      zwdw(:,:,1) = 0.e0
1751      DO jj = 1, jpj
1752         DO ji = 1, jpi
1753            DO jk = 2, jpk
1754               zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + & 
1755                              & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
1756            END DO
1757         END DO
1758      END DO
1759      !
1760      ! Conditionnal horizontal Laplacian diffusion:
1761      !---------------------------------------------
1762      IF ( ll_lapdiff_cond ) THEN
1763         !
1764         zwdw_b(:,:,1) = 0._wp
1765         DO jj = 1, jpj
1766            DO ji = 1, jpi
1767               DO jk=2,jpk
1768                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
1769                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
1770               END DO
1771            END DO
1772         END DO
1773         !
1774         DO jk = 2, jpkm1
1775            zwu(:,:) = 0._wp
1776            zwv(:,:) = 0._wp
1777            DO jj = 1, jpjm1
1778               DO ji = 1, fs_jpim1   ! vector opt.                 
1779                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
1780                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
1781                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
1782                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
1783               END DO
1784            END DO
1785            DO jj = 2, jpjm1
1786               DO ji = fs_2, fs_jpim1   ! vector opt.
1787                  ztmp1 = ( zwu(ji-1,jj  ) - zwu(ji,jj) &
1788                     &  +   zwv(ji  ,jj-1) - zwv(ji,jj) ) * r1_e1e2t(ji,jj)
1789                  zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e1e2t(ji,jj)), 0._wp)
1790                  ztmp = SIGN(zh2, ztmp1)
1791                  zeu2 = zhdiff * e1e2t(ji,jj)*e1e2t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj))
1792                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
1793               END DO
1794            END DO
1795         END DO         
1796         !
1797      ENDIF
1798
1799      ! Conditionnal horizontal Bilaplacian diffusion:
1800      !-----------------------------------------------
1801      IF ( ll_blpdiff_cond ) THEN
1802         !
1803         zwdw_b(:,:,1) = 0._wp
1804         DO jj = 1, jpj
1805            DO ji = 1, jpi
1806               DO jk = 2,jpkm1
1807                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + & 
1808                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
1809               END DO
1810            END DO
1811         END DO
1812         !
1813         DO jk = 2, jpkm1
1814            zwu(:,:) = 0._wp
1815            zwv(:,:) = 0._wp
1816            DO jj = 1, jpjm1
1817               DO ji = 1, fs_jpim1   ! vector opt.                 
1818                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
1819                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
1820                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
1821                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
1822               END DO
1823            END DO
1824            DO jj = 2, jpjm1
1825               DO ji = fs_2, fs_jpim1   ! vector opt.
1826                  zwdw_b(ji,jj,jk) = -( (zwu(ji-1,jj  ) - zwu(ji,jj)) &
1827                                &  +    (zwv(ji  ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj)
1828               END DO
1829            END DO
1830         END DO   
1831         !
1832         CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )     
1833         !
1834         DO jk = 2, jpkm1
1835            zwu(:,:) = 0._wp
1836            zwv(:,:) = 0._wp
1837            DO jj = 1, jpjm1
1838               DO ji = 1, fs_jpim1   ! vector opt.                 
1839                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) &
1840                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) )
1841                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) & 
1842                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) )
1843               END DO
1844            END DO
1845            DO jj = 2, jpjm1
1846               DO ji = fs_2, fs_jpim1   ! vector opt.
1847                  ztmp1 = ( (zwu(ji-1,jj  ) - zwv(ji,jj)) &
1848                     &  +   (zwv(ji  ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj)
1849                  zh2 = MAX(abs(ztmp1)-zhlim2*SQRT(r1_e1e2t(ji,jj))*r1_e1e2t(ji,jj), 0._wp)
1850                  ztmp = SIGN(zh2, ztmp1)
1851                  zeu2 = zhdiff2 * e1e2t(ji,jj)*e1e2t(ji,jj) / 16._wp
1852                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk)
1853               END DO
1854            END DO
1855         END DO   
1856         !
1857      ENDIF
1858
1859      ! Conditionnal vertical diffusion:
1860      !---------------------------------
1861      IF ( ll_zdiff_cond ) THEN
1862         DO jk = 2, jpkm1
1863            DO jj = 2, jpjm1
1864               DO ji = fs_2, fs_jpim1   ! vector opt.   
1865                  ztmp  = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) & 
1866                            -(tilde_e3t_b(ji,jj,jk  )+e3t_0(ji,jj,jk  ))*tmask(ji,jj,jk  ) ) 
1867                  ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) & 
1868                        &           +tilde_e3t_b(ji,jj,jk  ) + e3t_0(ji,jj,jk  ) )
1869                  zh2   = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp)
1870                  ztmp  = SIGN(zh2, ztmp)
1871                  IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0
1872                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk)
1873               END DO
1874            END DO
1875         END DO
1876      ENDIF
1877      !
1878      ! Check grid from the bottom to the surface
1879      !------------------------------------------
1880      IF ( ll_chk_bot2top ) THEN
1881         DO jj = 2, jpjm1
1882            DO ji = 2, jpim1
1883               jkbot = mbkt(ji,jj)                   
1884               DO jk = jkbot,2,-1
1885                  !
1886                  zh_0   = e3t_0(ji,jj,jk)
1887                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1))
1888                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
1889                  zh_min = MIN(zh_0/3._wp, dzmin_int)
1890!                  zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
1891                  !
1892                  ! Set maximum and minimum vertical excursions   
1893                  ztmph = hsm(ji,jj)
1894                  ztmpd = dsm(ji,jj)
1895                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)))/ztmpd)
1896!                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)
1897                  zh2 = MAX(zh2,0.001_wp) ! Extend tolerance a bit for stability reasons (to be explored)
1898                  zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 )
1899                  zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff)
1900                  zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 )
1901                  zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff)
1902                  !
1903                  ! New layer thickness:
1904                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
1905                  !
1906                  ! Ensure minimum layer thickness:                 
1907!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
1908!                  zh_new = cush(zh_new, zh_min)
1909                  zh_new = MAX(zh_new, zh_min)           
1910                  !
1911                  ! Final flux:
1912                  zdiff = (zh_new - zh_old) * tmask(ji,jj,jk)
1913                  !
1914                  ! Limit thickness change in 1 time step:
1915!                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
1916!                  zdiff = SIGN(ztmp, zh_new - zh_old)
1917                  zh_new = zdiff + zh_old
1918                  !
1919                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk+1) - zh_new       
1920               END DO   
1921            END DO
1922         END DO
1923      END IF   
1924      !
1925      ! Check grid from the surface to the bottom
1926      !------------------------------------------
1927      IF ( ll_chk_top2bot ) THEN     
1928         DO jj = 2, jpjm1
1929            DO ji = 2, jpim1
1930               jkbot = mbkt(ji,jj)   
1931               DO jk = 1, jkbot-1
1932                  !
1933                  zh_0   = e3t_0(ji,jj,jk)
1934                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1))
1935                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
1936                  zh_min = MIN(zh_0/3._wp, dzmin_int)
1937!                  zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
1938                  !
1939!                  zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf)
1940                  !
1941                  ! New layer thickness:
1942                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
1943                  !
1944                  ! Ensure minimum layer thickness:
1945!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
1946!                  zh_new = cush(zh_new, zh_min)
1947                  zh_new = MAX(zh_new, zh_min)   
1948                  !
1949                  ! Final flux:
1950                  zdiff = (zh_new -zh_old) * tmask(ji,jj,jk)
1951                  !
1952                  ! Limit flux:                 
1953!                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
1954!                  zdiff = SIGN(ztmp, zh_new - zh_old)
1955                  zh_new = zdiff + zh_old
1956                  !
1957                  zwdw(ji,jj,jk+1) = zwdw(ji,jj,jk) + zh_new
1958               END DO
1959               !               
1960            END DO
1961         END DO
1962      ENDIF
1963      !
1964      DO jj = 2, jpjm1
1965         DO ji = 2, jpim1
1966            DO jk = 1, jpkm1
1967               tilde_e3t_a(ji,jj,jk) =  (zwdw(ji,jj,jk+1)-zwdw(ji,jj,jk)-e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
1968            END DO
1969         END DO
1970      END DO
1971      !
1972      !
1973      IF( ln_timing )  CALL timing_stop('dom_vvl_regrid')
1974      !
1975   END SUBROUTINE dom_vvl_regrid
1976
1977   FUNCTION cush(hin, hmin)  RESULT(hout)
1978      !!----------------------------------------------------------------------
1979      !!                 ***  FUNCTION cush  ***
1980      !!
1981      !! ** Purpose :   
1982      !!
1983      !! ** Method  :
1984      !!
1985      !!----------------------------------------------------------------------
1986      IMPLICIT NONE
1987      REAL(wp), INTENT(in) ::  hin, hmin
1988      REAL(wp)             ::  hout, zx, zh_cri
1989      !!----------------------------------------------------------------------
1990      zh_cri = 3._wp * hmin
1991      !
1992      IF ( hin<=0._wp ) THEN
1993         hout = hmin
1994      !
1995      ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN
1996         zx = hin/zh_cri
1997         hout = hmin * (1._wp + zx + zx*zx)     
1998      !
1999      ELSEIF ( hin>zh_cri ) THEN
2000         hout = hin
2001      !
2002      ENDIF
2003      !
2004   END FUNCTION cush
2005
2006   FUNCTION cush_max(hin, hmax)  RESULT(hout)
2007      !!----------------------------------------------------------------------
2008      !!                 ***  FUNCTION cush  ***
2009      !!
2010      !! ** Purpose :   
2011      !!
2012      !! ** Method  :
2013      !!
2014      !!----------------------------------------------------------------------
2015      IMPLICIT NONE
2016      REAL(wp), INTENT(in) ::  hin, hmax
2017      REAL(wp)             ::  hout, hmin, zx, zh_cri
2018      !!----------------------------------------------------------------------
2019      hmin = 0.1_wp * hmax
2020      zh_cri = 3._wp * hmin
2021      !
2022      IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN
2023         zx = (hmax-hin)/zh_cri
2024         hout = hmax - hmin * (1._wp + zx + zx*zx)
2025         !
2026      ELSEIF ( hin>(hmax-zh_cri) ) THEN
2027         hout = hmax - hmin
2028         !
2029      ELSE
2030         hout = hin
2031         !
2032      ENDIF
2033      !
2034   END FUNCTION cush_max
2035
2036   SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin )
2037      !!----------------------------------------------------------------------
2038      !!                  ***  ROUTINE dom_vvl_adv_fct  ***
2039      !!
2040      !! **  Purpose :  Do thickness advection
2041      !!
2042      !! **  Method  :  FCT scheme to ensure positivity
2043      !!
2044      !! **  Action  : - Update pta thickness tendency and diffusive fluxes
2045      !!               - this is the total trend, hence it does include sea level motions
2046      !!----------------------------------------------------------------------
2047      !
2048      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
2049      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta    ! thickness baroclinic trend
2050      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input velocities
2051      !
2052      INTEGER  ::   ji, jj, jk, ib, ib_bdy               ! dummy loop indices 
2053      INTEGER  ::   ikbu, ikbv, ibot
2054      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
2055      REAL(wp) ::   zdi, zdj, zmin           !   -      -
2056      REAL(wp) ::   zfp_ui, zfp_vj           !   -      -
2057      REAL(wp) ::   zfm_ui, zfm_vj           !   -      -
2058      REAL(wp) ::   zfp_hi, zfp_hj           !   -      -
2059      REAL(wp) ::   zfm_hi, zfm_hj           !   -      -
2060      REAL(wp) ::   ztout,  ztin, zfac       !   -      -
2061      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwi
2062      !!----------------------------------------------------------------------
2063      !
2064      IF( ln_timing )  CALL timing_start('dom_vvl_adv_fct')
2065      !
2066      !
2067      ! 1. Initializations
2068      ! ------------------
2069      !
2070      IF( neuler == 0 .AND. kt == nit000 ) THEN
2071         z2dtt =  rdt
2072      ELSE
2073         z2dtt = 2.0_wp * rdt
2074      ENDIF
2075      !
2076      zwi(:,:,:) = 0.e0
2077      zwx(:,:,:) = 0.e0 
2078      zwy(:,:,:) = 0.e0
2079      !
2080      !
2081      ! 2. upstream advection with initial mass fluxes & intermediate update
2082      ! --------------------------------------------------------------------
2083      IF ( ll_shorizd ) THEN
2084         DO jk = 1, jpkm1
2085            DO jj = 1, jpjm1
2086               DO ji = 1, fs_jpim1   ! vector opt.
2087                  !
2088                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
2089                  zfp_hi = MIN(zfp_hi, e3t_b(ji  ,jj  ,jk))
2090                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
2091                  !
2092                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp)
2093                  zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj  ,jk))
2094                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 
2095                  !
2096                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
2097                  zfp_hj = MIN(zfp_hj, e3t_b(ji  ,jj  ,jk))
2098                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
2099                  !
2100                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp)
2101                  zfm_hj = MIN(zfm_hj, e3t_b(ji  ,jj+1,jk))
2102                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
2103                  !
2104                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
2105                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
2106                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
2107                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
2108                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2109                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2110               END DO
2111            END DO
2112         END DO
2113      ELSE
2114         DO jk = 1, jpkm1
2115            DO jj = 1, jpjm1
2116               DO ji = 1, fs_jpim1   ! vector opt.
2117                  !
2118                  zfp_hi = e3t_b(ji  ,jj  ,jk)
2119                  zfm_hi = e3t_b(ji+1,jj  ,jk)             
2120                  zfp_hj = e3t_b(ji  ,jj  ,jk)
2121                  zfm_hj = e3t_b(ji  ,jj+1,jk)
2122                  !
2123                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
2124                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
2125                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
2126                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
2127                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2128                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2129               END DO
2130            END DO
2131         END DO
2132      ENDIF
2133
2134      ! total advective trend
2135      DO jk = 1, jpkm1
2136         DO jj = 2, jpjm1
2137            DO ji = fs_2, fs_jpim1   ! vector opt.
2138               zbtr = r1_e1e2t(ji,jj)
2139               ! total intermediate advective trends
2140               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )  &
2141                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )  )
2142               !
2143               ! update and guess with monotonic sheme
2144               pta(ji,jj,jk) =   pta(ji,jj,jk) + ztra
2145               zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk)
2146            END DO
2147         END DO
2148      END DO
2149
2150      CALL lbc_lnk( zwi, 'T', 1. ) 
2151
2152      IF ( ln_bdy ) THEN
2153         DO ib_bdy=1, nb_bdy
2154            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
2155               ji = idx_bdy(ib_bdy)%nbi(ib,1)
2156               jj = idx_bdy(ib_bdy)%nbj(ib,1)
2157               DO jk = 1, jpkm1 
2158                  zwi(ji,jj,jk) = e3t_a(ji,jj,jk)
2159               END DO
2160            END DO
2161         END DO
2162      ENDIF
2163
2164!      IF ( ln_vvl_dbg ) THEN
2165!         zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 )
2166!         IF( lk_mpp )   CALL mpp_min( zmin )
2167!         IF( zmin < 0._wp) THEN
2168!            IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here')
2169!            IF(lwp) WRITE(numout,*) zmin
2170!         ENDIF
2171!      ENDIF
2172
2173      ! 3. antidiffusive flux : high order minus low order
2174      ! --------------------------------------------------
2175      ! antidiffusive flux on i and j
2176      DO jk = 1, jpkm1
2177         DO jj = 1, jpjm1
2178            DO ji = 1, fs_jpim1   ! vector opt.
2179               zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * e3u_n(ji,jj,jk) &
2180                                & - zwx(ji,jj,jk)) * umask(ji,jj,jk)
2181               zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * e3v_n(ji,jj,jk) & 
2182                                & - zwy(ji,jj,jk)) * vmask(ji,jj,jk)
2183               !
2184               ! Update advective fluxes
2185               un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk)
2186               vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk)
2187            END DO
2188         END DO
2189      END DO
2190     
2191      CALL lbc_lnk_multi( zwx, 'U', -1., zwy, 'V', -1. )     !* local domain boundaries
2192
2193      ! 4. monotonicity algorithm
2194      ! -------------------------
2195      CALL nonosc_2d( e3t_b(:,:,:), zwx, zwy, zwi, z2dtt )
2196
2197      ! 5. final trend with corrected fluxes
2198      ! ------------------------------------
2199      !
2200      ! Update advective fluxes
2201      un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:)
2202      vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:)
2203      !
2204      DO jk = 1, jpkm1
2205         DO jj = 2, jpjm1
2206            DO ji = fs_2, fs_jpim1   ! vector opt. 
2207               !
2208               zbtr = r1_e1e2t(ji,jj)
2209               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
2210                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
2211               ! add them to the general tracer trends
2212               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
2213               zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * pta(ji,jj,jk)* bdytmask(ji,jj) ) * tmask(ji,jj,jk)
2214            END DO
2215         END DO
2216      END DO
2217
2218      IF ( ln_vvl_dbg ) THEN
2219         zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 ) 
2220         IF( lk_mpp )   CALL mpp_min( zmin )
2221         IF( zmin < 0._wp) THEN
2222            IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here')
2223            IF(lwp) WRITE(numout,*) zmin
2224         ENDIF
2225      ENDIF
2226      !
2227      IF( ln_timing )  CALL timing_stop('dom_vvl_adv_fct')
2228      !
2229   END SUBROUTINE dom_vvl_adv_fct
2230
2231   SUBROUTINE dom_vvl_ups_cor( kt, pta, uin, vin ) 
2232      !!----------------------------------------------------------------------
2233      !!                  ***  ROUTINE dom_vvl_adv_fct  ***
2234      !!
2235      !! **  Purpose :  Correct for addionnal barotropic fluxes
2236      !!                in the upstream direction
2237      !!
2238      !! **  Method  :   
2239      !!
2240      !! **  Action  : - Update diffusive fluxes uin, vin
2241      !!               - Remove divergence from thickness tendency
2242      !!----------------------------------------------------------------------
2243      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
2244      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta       ! thickness baroclinic trend
2245      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input fluxes
2246      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
2247      INTEGER  ::   ikbu, ikbv, ibot
2248      REAL(wp) ::   zbtr, ztra               ! local scalar
2249      REAL(wp) ::   zdi, zdj                 !   -      -
2250      REAL(wp) ::   zfp_hi, zfp_hj           !   -      -
2251      REAL(wp) ::   zfm_hi, zfm_hj           !   -      -
2252      REAL(wp) ::   zfp_ui, zfp_vj           !   -      -
2253      REAL(wp) ::   zfm_ui, zfm_vj           !   -      -
2254      REAL(wp), DIMENSION(jpi,jpj) :: zbu, zbv, zhu_b, zhv_b
2255      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy
2256      !!----------------------------------------------------------------------
2257      !
2258      IF( ln_timing )  CALL timing_start('dom_vvl_ups_cor')
2259      !
2260      ! Compute barotropic flux difference:
2261      zbu(:,:) = 0.e0
2262      zbv(:,:) = 0.e0
2263      DO jj = 1, jpj
2264         DO ji = 1, jpi   ! vector opt.
2265            DO jk = 1, jpkm1
2266               zbu(ji,jj) = zbu(ji,jj) - uin(ji,jj,jk) * umask(ji,jj,jk)
2267               zbv(ji,jj) = zbv(ji,jj) - vin(ji,jj,jk) * vmask(ji,jj,jk)
2268            END DO
2269         END DO
2270      ENDDO
2271 
2272      ! Compute upstream depths:
2273      zhu_b(:,:) = 0.e0
2274      zhv_b(:,:) = 0.e0
2275
2276      IF ( ll_shorizd ) THEN
2277         ! Correct bottom value
2278         ! considering "shelf horizon depth"
2279         DO jj = 1, jpjm1
2280            DO ji = 1, jpim1   ! vector opt.
2281               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
2282               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
2283               DO jk=1, jpkm1
2284                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
2285                  zfp_hi = MIN(zfp_hi, e3t_b(ji  ,jj  ,jk))
2286                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
2287                  !
2288                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp)
2289                  zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj  ,jk))
2290                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
2291                  !
2292                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
2293                  zfp_hj = MIN(zfp_hj, e3t_b(ji  ,jj  ,jk))
2294                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
2295                  !
2296                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp)
2297                  zfm_hj = MIN(zfm_hj, e3t_b(ji  ,jj+1,jk))
2298                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
2299                  !
2300                  zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi &
2301                               &                 + (1._wp-zdi) * zfm_hi &
2302                               &                ) * umask(ji,jj,jk)
2303                  zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj &
2304                               &                 + (1._wp-zdj) * zfm_hj &
2305                               &                ) * vmask(ji,jj,jk)
2306               END DO
2307            END DO
2308         END DO
2309      ELSE
2310         DO jj = 1, jpjm1
2311            DO ji = 1, jpim1   ! vector opt.
2312               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 
2313               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
2314               DO jk = 1, jpkm1
2315                  zfp_hi = e3t_b(ji  ,jj  ,jk)
2316                  zfm_hi = e3t_b(ji+1,jj  ,jk)
2317                  zfp_hj = e3t_b(ji  ,jj  ,jk)
2318                  zfm_hj = e3t_b(ji  ,jj+1,jk)
2319                  !
2320                  zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi &
2321                               &                  + (1._wp-zdi) * zfm_hi &
2322                               &                ) * umask(ji,jj,jk)
2323                  !
2324                  zhv_b(ji,jj) = zhv_b(ji,jj) + (          zdj  * zfp_hj &
2325                               &                  + (1._wp-zdj) * zfm_hj &
2326                               &                ) * vmask(ji,jj,jk)
2327               END DO
2328            END DO
2329         END DO
2330      ENDIF
2331
2332      CALL lbc_lnk_multi( zhu_b(:,:), 'U', 1., zhv_b(:,:), 'V', 1. )     !* local domain boundaries
2333
2334      ! Corrective barotropic velocity (times hor. scale factor)
2335      zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1))
2336      zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1))
2337     
2338      ! Set corrective fluxes in upstream direction:
2339      !
2340      zwx(:,:,:) = 0.e0
2341      zwy(:,:,:) = 0.e0
2342
2343      IF ( ll_shorizd ) THEN
2344         DO jj = 1, jpjm1
2345            DO ji = 1, fs_jpim1   ! vector opt.
2346               ! upstream scheme
2347               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
2348               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
2349               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
2350               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
2351               DO jk = 1, jpkm1
2352                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
2353                  zfp_hi = MIN(e3t_b(ji  ,jj  ,jk), zfp_hi)
2354                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 
2355                  !
2356                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp)
2357                  zfm_hi = MIN(e3t_b(ji+1,jj  ,jk), zfm_hi)
2358                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 
2359                  !
2360                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp)
2361                  zfp_hj = MIN(e3t_b(ji  ,jj  ,jk), zfp_hj) 
2362                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 
2363                  !
2364                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp)
2365                  zfm_hj = MIN(e3t_b(ji  ,jj+1,jk), zfm_hj)
2366                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 
2367                  !
2368                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2369                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2370               END DO
2371            END DO
2372         END DO
2373      ELSE
2374         DO jj = 1, jpjm1
2375            DO ji = 1, fs_jpim1   ! vector opt.
2376               ! upstream scheme
2377               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
2378               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
2379               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
2380               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
2381               DO jk = 1, jpkm1
2382                  zfp_hi = e3t_b(ji  ,jj  ,jk)
2383                  zfm_hi = e3t_b(ji+1,jj  ,jk)
2384                  zfp_hj = e3t_b(ji  ,jj  ,jk)
2385                  zfm_hj = e3t_b(ji  ,jj+1,jk)
2386                  !
2387                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2388                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2389               END DO
2390            END DO
2391         END DO
2392      ENDIF
2393
2394      CALL lbc_lnk_multi( zwx, 'U', -1., zwy, 'V', -1. )     !* local domain boundaries
2395
2396      uin(:,:,:) = uin(:,:,:) + zwx(:,:,:)
2397      vin(:,:,:) = vin(:,:,:) + zwy(:,:,:)
2398      !
2399      ! Update trend with corrective fluxes:
2400      DO jk = 1, jpkm1
2401         DO jj = 2, jpjm1
2402            DO ji = fs_2, fs_jpim1   ! vector opt. 
2403               !
2404               zbtr = r1_e1e2t(ji,jj)
2405               ! total advective trends
2406               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
2407                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
2408               ! add them to the general tracer trends
2409               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
2410            END DO
2411         END DO
2412      END DO
2413      !
2414      IF( ln_timing )  CALL timing_stop('dom_vvl_ups_cor')
2415      !
2416   END SUBROUTINE dom_vvl_ups_cor
2417
2418   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt )
2419      !!---------------------------------------------------------------------
2420      !!                    ***  ROUTINE nonosc_2d  ***
2421      !!     
2422      !! **  Purpose :   compute monotonic thickness fluxes from the upstream
2423      !!       scheme and the before field by a nonoscillatory algorithm
2424      !!
2425      !! **  Method  :   ... ???
2426      !!       warning : pbef and paft must be masked, but the boundaries
2427      !!       conditions on the fluxes are not necessary zalezak (1979)
2428      !!       drange (1995) multi-dimensional forward-in-time and upstream-
2429      !!       in-space based differencing for fluid
2430      !!----------------------------------------------------------------------
2431      !
2432      !!----------------------------------------------------------------------
2433      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
2434      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
2435      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb        ! monotonic fluxes in the 3 directions
2436      !
2437      INTEGER ::   ji, jj, jk   ! dummy loop indices
2438      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
2439      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
2440      REAL(wp) ::   zupip1, zupim1, zupjp1, zupjm1, zupb, zupa
2441      REAL(wp) ::   zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa
2442      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
2443      !!----------------------------------------------------------------------
2444      !
2445      IF( ln_timing )  CALL timing_start('nonosc2')
2446      !
2447      zbig  = 1.e+40_wp
2448      zrtrn = 1.e-15_wp
2449      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp
2450
2451
2452      ! Search local extrema
2453      ! --------------------
2454      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
2455      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   &
2456         &        paft * tmask - zbig * ( 1.e0 - tmask )  )
2457      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   &
2458         &        paft * tmask + zbig * ( 1.e0 - tmask )  )
2459
2460      DO jk = 1, jpkm1
2461         z2dtt = p2dt
2462         DO jj = 2, jpjm1
2463            DO ji = fs_2, fs_jpim1   ! vector opt.
2464
2465               ! search maximum in neighbourhood
2466               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
2467                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
2468                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ))
2469
2470               ! search minimum in neighbourhood
2471               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
2472                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
2473                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ))
2474
2475               ! positive part of the flux
2476               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
2477                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   
2478
2479               ! negative part of the flux
2480               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
2481                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )
2482
2483               ! up & down beta terms
2484               zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt
2485               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
2486               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
2487            END DO
2488         END DO
2489      END DO
2490      CALL lbc_lnk_multi( zbetup, 'T', 1. , zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
2491
2492      ! 3. monotonic flux in the i & j direction (paa & pbb)
2493      ! ----------------------------------------
2494      DO jk = 1, jpkm1
2495         DO jj = 2, jpjm1
2496            DO ji = fs_2, fs_jpim1   ! vector opt.
2497               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
2498               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
2499               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
2500               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
2501
2502               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
2503               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
2504               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
2505               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
2506            END DO
2507         END DO
2508      END DO
2509      CALL lbc_lnk_multi( paa, 'U', -1., pbb, 'V', -1. )     !* local domain boundaries
2510      !
2511      IF( ln_timing )  CALL timing_stop('nonosc2')
2512      !
2513   END SUBROUTINE nonosc_2d
2514
2515   !!======================================================================
2516END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.