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

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

ztilde update, #2126: Add diagnostics

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