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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10116 – NEMO

Changeset 10116


Ignore:
Timestamp:
2018-09-12T17:39:01+02:00 (6 years ago)
Author:
jchanut
Message:

ztilde update (3): #2126
Rewrite thickness advection, add regriding, biharmonic interface diffusion

Location:
NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA/diawri.F90

    r9652 r10116  
    125125      ! Output of initial vertical scale factor 
    126126      CALL iom_put("e3t_0", e3t_0(:,:,:) ) 
    127       CALL iom_put("e3u_0", e3t_0(:,:,:) ) 
    128       CALL iom_put("e3v_0", e3t_0(:,:,:) ) 
     127      CALL iom_put("e3u_0", e3u_0(:,:,:) ) 
     128      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
    129129      ! 
    130130      CALL iom_put( "e3t" , e3t_n(:,:,:) ) 

  • NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90

    r10060 r10116  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 
    99   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
     10   !!            4.0  !  2018-09  (J. Chanut) improve z_tilde robustness 
    1011   !!---------------------------------------------------------------------- 
    1112 
    … …  
    3132   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3233   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 
    3337 
    3438   IMPLICIT NONE 
    … …  
    4448   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate 
    4549   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate 
    46    LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate 
    47    LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate 
    48    LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer 
     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 
    4962   !                                                       ! conservation: not used yet 
    50    REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient 
    51    REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days] 
    52    REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days] 
    53    REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation 
    54    LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints 
     63   REAL(wp)         :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian) 
     64   REAL(wp)         :: rn_ahe3_blp               ! thickness diffusion coefficient (Bilaplacian) 
     65   REAL(wp)         :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days] 
     66   REAL(wp)         :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days] 
     67   REAL(wp)         :: rn_day_ramp               ! Duration of linear ramp  [days] 
     68   REAL(wp)         :: hsmall=0.01_wp            ! small thickness [m] 
    5569 
    5670   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport 
    57    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence 
     71   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_lf, vn_lf, hdivn_lf      ! low frequency fluxes and divergence 
    5872   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors 
    5973   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors 
    60    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors 
    61    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
     74   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: tildemask                   ! mask tilde tendency 
     75   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! restoring period for scale factors 
     76   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! restoring period for low freq. divergence 
     77   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: hsm, dsm                    !  
     78   INTEGER         , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: i_int_bot 
    6279 
    6380   !! * Substitutions 
    … …  
    7895         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   & 
    7996            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   & 
    80             &      STAT = dom_vvl_alloc        ) 
     97            &      tildemask(jpi,jpj) , hsm(jpi,jpj) , dsm(jpi,jpj) , i_int_bot(jpi,jpj), STAT = dom_vvl_alloc ) 
    8198         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
    8299         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 
    … …  
    85102      ENDIF 
    86103      IF( ln_vvl_ztilde ) THEN 
    87          ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 
     104         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdivn_lf(jpi,jpj,jpk)             ,    & 
     105            &      un_lf(jpi,jpj,jpk), vn_lf(jpi,jpj,jpk)      , STAT= dom_vvl_alloc ) 
    88106         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
    89107         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 
    … …  
    117135      INTEGER ::   ji, jj, jk 
    118136      INTEGER ::   ii0, ii1, ij0, ij1 
    119       REAL(wp)::   zcoef 
     137      REAL(wp)::   zcoef, zwgt, ztmp, zhmin, zhmax 
    120138      !!---------------------------------------------------------------------- 
    121139      ! 
    … …  
    129147      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    130148      ! 
    131       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
     149      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdivn_lf 
    132150      CALL dom_vvl_rst( nit000, 'READ' ) 
    133151      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    … …  
    199217 
    200218      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies) 
     219      tildemask(:,:) = 1._wp 
     220 
    201221      IF( ln_vvl_ztilde ) THEN 
    202222!!gm : idea: add here a READ in a file of custumized restoring frequency 
    203          !                                   ! Values in days provided via the namelist 
    204          !                                   ! use rsmall to avoid possible division by zero errors with faulty settings 
    205          frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp ) 
    206          frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
    207          ! 
    208          IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile 
    209             frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings 
    210             frq_rst_hdv(:,:) = 1._wp / rdt 
     223         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 
     224         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp ) 
     225         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
     226         ! 
     227         IF( ln_vvl_ztilde_as_zstar ) THEN 
     228            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 
     229            frq_rst_e3t(:,:) = 0.0_wp  
     230            frq_rst_hdv(:,:) = 1.0_wp / rdt 
     231            tildemask(:,:) = 0._wp 
    211232         ENDIF 
    212          IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
     233         
     234         IF ( ln_vvl_zstar_at_eqtor ) THEN 
    213235            DO jj = 1, jpj 
    214236               DO ji = 1, jpi 
    215 !!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
     237!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default   
    216238                  IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    217239                     ! values outside the equatorial band and transition zone (ztilde) 
    218240                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    219                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
    220                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
     241!                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     242               
     243                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 
    221244                     ! values inside the equatorial band (ztilde as zstar) 
    222245                     frq_rst_e3t(ji,jj) =  0.0_wp 
    223                      frq_rst_hdv(ji,jj) =  1.0_wp / rdt 
    224                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    225                      !                                      ! (linearly transition from z-tilde to z-star) 
     246!                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt 
     247                     tildemask(ji,jj) = 0._wp 
     248                  ELSE 
     249                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values) 
    226250                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    227251                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    228252                        &                                          * 180._wp / 3.5_wp ) ) 
    229                      frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                & 
    230                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   & 
    231                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    232                         &                                          * 180._wp / 3.5_wp ) ) 
     253!                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                & 
     254!                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   & 
     255!                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     256!                        &                                          * 180._wp / 3.5_wp ) ) 
     257                     tildemask(ji,jj) = 0.5_wp * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     258                        &                                                 * 180._wp / 3.5_wp ) ) 
    233259                  ENDIF 
    234260               END DO 
    235261            END DO 
    236             IF( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
    237                ii0 = 103   ;   ii1 = 111        
    238                ij0 = 128   ;   ij1 = 135   ;    
    239                frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp 
    240                frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt 
    241             ENDIF 
    242262         ENDIF 
     263         ! 
     264         IF ( ln_vvl_zstar_on_shelf ) THEN 
     265            zhmin =  50._wp 
     266            zhmax = 100._wp 
     267            DO jj = 1, jpj 
     268               DO ji = 1, jpi 
     269                  zwgt = 1._wp 
     270                  IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN 
     271                     zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin) 
     272                  ELSEIF ( ht_0(ji,jj)<=zhmin) THEN 
     273                     zwgt = 0._wp 
     274                  ENDIF 
     275                  frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt) 
     276                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt) 
     277               END DO 
     278            END DO 
     279         ENDIF 
     280         ! 
     281         ztmp = MAXVAL( frq_rst_hdv(:,:) ) 
     282         IF( lk_mpp )   CALL mpp_max( ztmp )                 ! max over the global domain 
     283         ! 
     284         IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' ) 
     285         ! 
    243286      ENDIF 
    244287      ! 
    … …  
    256299         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    257300            !                                        ! ------------ ! 
    258             CALL iom_set_rstw_var_active('hdiv_lf') 
     301            CALL iom_set_rstw_var_active('un_lf') 
     302            CALL iom_set_rstw_var_active('vn_lf') 
     303            CALL iom_set_rstw_var_active('hdivn_lf') 
    259304         ENDIF 
    260305         ! 
    … …  
    279324      !!                               to the "baroclinic" level thickness. 
    280325      !! 
    281       !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case 
     326      !! ** Action  :  - hdivn_lf   : restoring towards full baroclinic divergence in z_tilde case 
    282327      !!               - tilde_e3t_a: after increment of vertical scale factor  
    283328      !!                              in z_tilde case 
    … …  
    289334      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence 
    290335      ! 
     336      LOGICAL                ::   ll_do_bclinic         ! local logical 
    291337      INTEGER                ::   ji, jj, jk            ! dummy loop indices 
     338      INTEGER                ::   ib, ib_bdy, ip, jp    !   "     "     " 
    292339      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers 
    293       REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars 
    294       LOGICAL                ::   ll_do_bclinic         ! local logical 
     340      INTEGER                ::   ncall 
     341      REAL(wp)               ::   z2dt  , z_tmin, z_tmax! local scalars                 
     342      REAL(wp)               ::   zalpha, zwgt          ! temporary scalars 
     343      REAL(wp)               ::   zdu, zdv, zramp 
     344      REAL(wp)               ::   ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv 
    295345      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv 
    296       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t 
     346      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t, ztu, ztv 
    297347      !!---------------------------------------------------------------------- 
    298348      ! 
    … …  
    308358 
    309359      ll_do_bclinic = .TRUE. 
     360      ncall = 1 
     361 
    310362      IF( PRESENT(kcall) ) THEN 
    311          IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE. 
     363         ! comment line below if tilda coordinate has to be computed at each call 
     364         ! This is not working yet 
     365!         IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE. 
     366         ncall = kcall   
     367      ENDIF 
     368 
     369      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     370         z2dt =  rdt 
     371      ELSE 
     372         z2dt = 2.0_wp * rdt 
    312373      ENDIF 
    313374 
    314375      ! ******************************* ! 
    315       ! After acale factors at t-points ! 
     376      ! After scale factors at t-points ! 
    316377      ! ******************************* ! 
    317378      !                                                ! --------------------------------------------- ! 
    … …  
    325386      END DO 
    326387      ! 
    327       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    328          !                                                            ! ------baroclinic part------ ! 
    329          ! I - initialization 
    330          ! ================== 
    331  
    332          ! 1 - barotropic divergence 
    333          ! ------------------------- 
    334          zhdiv(:,:) = 0._wp 
    335          zht(:,:)   = 0._wp 
     388      IF( (ln_vvl_ztilde.OR.ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     389         !                                                               ! ------baroclinic part------ ! 
     390         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms 
     391         un_td(:,:,:) = 0.0_wp        ! Transport corrections 
     392         vn_td(:,:,:) = 0.0_wp 
     393 
     394         zhdiv(:,:) = 0. 
    336395         DO jk = 1, jpkm1 
    337396            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
    338             zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    339          END DO 
    340          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
    341  
    342          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only) 
    343          ! -------------------------------------------------- 
    344          IF( ln_vvl_ztilde ) THEN 
    345             IF( kt > nit000 ) THEN 
     397         END DO 
     398         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )  
     399 
     400         ze3t(:,:,:) = 0._wp 
     401         IF( ln_rnf ) THEN 
     402            CALL sbc_rnf_div( ze3t )          ! runoffs  
     403            DO jk=1,jpkm1 
     404               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ze3t(:,:,jk) * e3t_n(:,:,jk) 
     405            END DO    
     406         ENDIF  
     407 
     408         ! Thickness advection: 
     409         ! -------------------- 
     410         ! Set advection velocities and source term 
     411         IF ( ln_vvl_ztilde ) THEN 
     412            IF ( ncall==1 ) THEN 
     413               zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
    346414               DO jk = 1, jpkm1 
    347                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
    348                      &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
    349                END DO 
     415                  ztu(:,:,jk) = un(:,:,jk) * e3u_n(:,:,jk) * e2u(:,:) 
     416                  ztv(:,:,jk) = vn(:,:,jk) * e3v_n(:,:,jk) * e1v(:,:) 
     417                  ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * zhdiv(:,:) 
     418               END DO 
     419               ! 
     420               un_lf(:,:,:)  =    un_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:) 
     421               vn_lf(:,:,:)  =    vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:) 
     422            hdivn_lf(:,:,:)  = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:) 
     423! 2nd order filtering: 
     424!               un_lf2(:,:,:)  =    un_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:) 
     425!               vn_lf2(:,:,:)  =    vn_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:) 
     426!            hdivn_lf2(:,:,:)  = hdivn_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:) 
     427!               un_lf(:,:,:)   =     un_lf(:,:,:) * (1._wp - zalpha) + zalpha * un_lf2(:,:,:) 
     428!               vn_lf(:,:,:)   =     vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * vn_lf2(:,:,:) 
     429!            hdivn_lf(:,:,:)   =  hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * hdivn_lf2(:,:,:) 
    350430            ENDIF 
     431            ! 
     432            DO jk = 1, jpkm1 
     433               ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/e3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk) 
     434               ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/e3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk) 
     435               tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
     436            END DO 
     437 
     438            ! 
     439         ELSEIF ( ln_vvl_layer ) THEN 
     440            ! 
     441            DO jk = 1, jpkm1 
     442               ztu(:,:,jk) = un(:,:,jk) 
     443               ztv(:,:,jk) = vn(:,:,jk) 
     444            END DO 
     445            ! 
    351446         ENDIF 
    352  
    353          ! II - after z_tilde increments of vertical scale factors 
    354          ! ======================================================= 
    355          tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms 
    356  
    357          ! 1 - High frequency divergence term 
    358          ! ---------------------------------- 
    359          IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
     447         ! 
     448         ! Block fluxes through small layers: 
     449         DO jk=1,jpkm1 
     450            DO ji = 1, jpi 
     451               DO jj= 1, jpj 
     452                  zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_n(ji,jj,jk) - hsmall) ) 
     453                  un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * e3u_n(ji,jj,jk) * e2u(ji,jj) 
     454                  ztu(ji,jj,jk) = zmsku * ztu(ji,jj,jk) 
     455                  IF ( ln_vvl_ztilde ) un_lf(ji,jj,jk) = zmsku * un_lf(ji,jj,jk) 
     456                  ! 
     457                  zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_n(ji,jj,jk) - hsmall) ) 
     458                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * e3v_n(ji,jj,jk) * e1v(ji,jj) 
     459                  ztv(ji,jj,jk) = zmskv * ztv(ji,jj,jk) 
     460                  IF ( ln_vvl_ztilde ) vn_lf(ji,jj,jk) = zmskv * vn_lf(ji,jj,jk) 
     461               END DO 
     462            END DO 
     463         END DO       
     464         ! 
     465         ! Do advection 
     466         IF     (ln_vvl_adv_fct) THEN 
     467            CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv ) 
     468            ! 
     469         ELSEIF (ln_vvl_adv_cn2) THEN 
    360470            DO jk = 1, jpkm1 
    361                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    362             END DO 
    363          ELSE                         ! layer case 
     471               DO jj = 2, jpjm1 
     472                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     473                     tilde_e3t_a(ji,jj,jk) =   tilde_e3t_a(ji,jj,jk) & 
     474                     & -(  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)       & 
     475                     &   + 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)  )    & 
     476                     & / ( e1t(ji,jj) * e2t(ji,jj) ) 
     477                  END DO 
     478               END DO 
     479            END DO 
     480         ENDIF 
     481         !  
     482         ! Thickness anomaly diffusion: 
     483         ! ----------------------------- 
     484         ztu(:,:,:) = 0.0_wp 
     485         ztv(:,:,:) = 0.0_wp 
     486 
     487         ze3t(:,:,1) = 0.e0 
     488         DO jj = 1, jpj 
     489            DO ji = 1, jpi 
     490               DO jk = 2, jpk 
     491                  ze3t(ji,jj,jk) =  ze3t(ji,jj,jk-1) + tilde_e3t_b(ji,jj,jk-1) * tmask(ji,jj,jk-1)  
     492               END DO 
     493            END DO 
     494         END DO 
     495 
     496         IF ( ln_vvl_blp ) THEN  ! Bilaplacian 
    364497            DO jk = 1, jpkm1 
    365                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
     498               DO jj = 1, jpjm1                 ! First derivative (gradient) 
     499                  DO ji = 1, fs_jpim1   ! vector opt.                   
     500!                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * e2_e1u(ji,jj) & 
     501!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     502!                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) &  
     503!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     504                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * e2_e1u(ji,jj) & 
     505                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) ) 
     506                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) &  
     507                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) ) 
     508                  END DO 
     509               END DO 
     510 
     511               DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
     512                  DO ji = fs_2, fs_jpim1 ! vector opt. 
     513                     zht(ji,jj) =  rn_ahe3_blp * r1_e1e2t(ji,jj) * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
     514                        &                                            + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   ) 
     515                  END DO 
     516               END DO 
     517 
     518               ! Open boundary conditions: 
     519               IF ( ln_bdy ) THEN 
     520                  DO ib_bdy=1, nb_bdy 
     521                     DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 
     522                        ji = idx_bdy(ib_bdy)%nbi(ib,1) 
     523                        jj = idx_bdy(ib_bdy)%nbj(ib,1) 
     524                        zht(ji,jj) = 0._wp 
     525                     END DO  
     526                  END DO 
     527               END IF 
     528 
     529               CALL lbc_lnk( zht, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn) 
     530 
     531               DO jj = 1, jpjm1                 ! third derivative (gradient) 
     532                  DO ji = 1, fs_jpim1   ! vector opt. 
     533                     ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) ) 
     534                     ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) * ( zht(ji  ,jj+1) - zht(ji,jj) ) 
     535                  END DO 
     536               END DO 
    366537            END DO 
    367538         ENDIF 
    368539 
    369          ! 2 - Restoring term (z-tilde case only) 
    370          ! ------------------ 
    371          IF( ln_vvl_ztilde ) THEN 
    372             DO jk = 1, jpk 
    373                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
    374             END DO 
    375          ENDIF 
    376  
    377          ! 3 - Thickness diffusion term 
    378          ! ---------------------------- 
    379          zwu(:,:) = 0._wp 
    380          zwv(:,:) = 0._wp 
    381          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    382             DO jj = 1, jpjm1 
    383                DO ji = 1, fs_jpim1   ! vector opt. 
    384                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    385                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    386                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    387                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    388                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    389                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    390                END DO 
    391             END DO 
    392          END DO 
    393          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    394             DO ji = 1, jpi 
    395                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    396                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    397             END DO 
    398          END DO 
    399          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
     540         IF ( ln_vvl_lap ) THEN    ! Laplacian 
     541            DO jk = 1, jpkm1                    ! First derivative (gradient) 
     542               DO jj = 1, jpjm1 
     543                  DO ji = 1, fs_jpim1   ! vector opt.                   
     544!                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) & 
     545!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     546!                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) &  
     547!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     548                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) & 
     549                         &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) ) 
     550                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) &  
     551                         &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) ) 
     552                     ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu 
     553                     ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv 
     554                  END DO 
     555               END DO 
     556            END DO 
     557         ENDIF  
     558 
     559         ! divergence of diffusive fluxes 
     560         DO jk = 1, jpkm1 
    400561            DO jj = 2, jpjm1 
    401562               DO ji = fs_2, fs_jpim1   ! vector opt. 
    402                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    403                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
     563!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk) - ztu(ji,jj,jk)    & 
     564!                     &                                          +     ztv(ji  ,jj-1,jk) - ztv(ji,jj,jk)    & 
     565!                     &                                            ) * r1_e1e2t(ji,jj) 
     566                  un_td(ji,jj,jk) = un_td(ji,jj,jk) + ztu(ji,jj,jk+1) - ztu(ji,jj,jk  )  
     567                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) + ztv(ji,jj,jk+1) - ztv(ji,jj,jk  )  
     568                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk+1) - ztu(ji,jj,jk+1)    & 
     569                     &                                               +ztv(ji  ,jj-1,jk+1) - ztv(ji,jj,jk+1)    & 
     570                     &                                               -ztu(ji-1,jj  ,jk  ) + ztu(ji,jj,jk  )    & 
     571                     &                                               -ztv(ji  ,jj-1,jk  ) + ztv(ji,jj,jk  )    & 
    404572                     &                                            ) * r1_e1e2t(ji,jj) 
    405573               END DO 
    406574            END DO 
    407575         END DO 
    408          !                       ! d - thickness diffusion transport: boundary conditions 
    409          !                             (stored for tracer advction and continuity equation) 
    410          CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp) 
    411  
    412          ! 4 - Time stepping of baroclinic scale factors 
    413          ! --------------------------------------------- 
     576 
     577  
     578!         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:) 
     579!         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:) 
     580 
     581         CALL lbc_lnk( un_td , 'U' , -1.) 
     582         CALL lbc_lnk( vn_td , 'V' , -1.)  
     583         ! 
     584         CALL dom_vvl_ups_cor( kt, tilde_e3t_a, un_td, vn_td ) 
     585 
     586!         IF ( ln_vvl_ztilde ) THEN 
     587!            ztu(:,:,:) = -un_lf(:,:,:) 
     588!            ztv(:,:,:) = -vn_lf(:,:,:) 
     589!            CALL dom_vvl_ups_cor( kt, tilde_e3t_a, ztu, ztv ) 
     590!         ENDIF 
     591         ! 
     592         ! Remove "external thickness" tendency: 
     593         DO jk = 1, jpkm1 
     594            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   e3t_n(:,:,jk) * zhdiv(:,:) 
     595         END DO  
     596                    
    414597         ! Leapfrog time stepping 
    415598         ! ~~~~~~~~~~~~~~~~~~~~~~ 
    416          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    417             z2dt =  rdt 
    418          ELSE 
    419             z2dt = 2.0_wp * rdt 
     599         zramp = 1._wp 
     600         IF ((.NOT.ln_rstart).AND.ln_vvl_ramp) zramp = MIN(MAX( ((kt-nit000)*rdt)/(rn_day_ramp*rday),0._wp),1._wp) 
     601 
     602         DO jk=1,jpkm1 
     603            tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) & 
     604                               & * tildemask(:,:) * zramp 
     605         END DO 
     606 
     607         ! Ensure layer separation: 
     608         ! ~~~~~~~~~~~~~~~~~~~~~~~~ 
     609         IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt ) 
     610   
     611         ! Boundary conditions: 
     612         ! ~~~~~~~~~~~~~~~~~~~~ 
     613         IF ( ln_bdy ) THEN 
     614            DO ib_bdy = 1, nb_bdy 
     615               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 
     616!!               DO ib = 1, idx_bdy(ib_bdy)%nblen(1) 
     617                  ji = idx_bdy(ib_bdy)%nbi(ib,1) 
     618                  jj = idx_bdy(ib_bdy)%nbj(ib,1) 
     619                  zwgt = idx_bdy(ib_bdy)%nbw(ib,1) 
     620                  ip = bdytmask(ji+1,jj  ) - bdytmask(ji-1,jj  ) 
     621                  jp = bdytmask(ji  ,jj+1) - bdytmask(ji  ,jj-1) 
     622                  DO jk = 1, jpkm1   
     623                     tilde_e3t_a(ji,jj,jk) = 0.e0 
     624!!                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt) 
     625!!                     tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk) 
     626                  END DO 
     627               END DO  
     628            END DO 
    420629         ENDIF 
    421          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    422          tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
    423  
    424          ! Maximum deformation control 
    425          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    426          ze3t(:,:,jpk) = 0._wp 
     630 
     631         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )         
     632         
     633      ENDIF 
     634 
     635      IF( ln_vvl_ztilde.AND.( ncall==1))  THEN 
     636         zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
     637         ! 
     638         ! divergence of diffusive fluxes 
    427639         DO jk = 1, jpkm1 
    428             ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 
    429          END DO 
    430          z_tmax = MAXVAL( ze3t(:,:,:) ) 
    431          IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain 
    432          z_tmin = MINVAL( ze3t(:,:,:) ) 
     640            DO jj = 2, jpjm1 
     641               DO ji = fs_2, fs_jpim1   ! vector opt. 
     642                  ze3t(ji,jj,jk) = (  un_td(ji,jj,jk) - un_td(ji-1,jj  ,jk) & 
     643                                 &  + vn_td(ji,jj,jk) - vn_td(ji  ,jj-1,jk) & 
     644                                 & ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     645               END DO 
     646            END DO 
     647         END DO 
     648         CALL lbc_lnk( ze3t(:,:,:), 'T', 1. ) 
     649         hdivn_lf(:,:,:)  =  hdivn_lf(:,:,:)  + zalpha * ze3t(:,:,:)  
     650! 2nd order filtering: 
     651!         hdivn_lf2(:,:,:)  = hdivn_lf2(:,:,:) + zalpha * ze3t(:,:,:) 
     652!         hdivn_lf(:,:,:)  =  hdivn_lf(:,:,:)  + zalpha * zalpha * ze3t(:,:,:) 
     653      ENDIF 
     654 
     655      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
     656      !                                             ! ---baroclinic part--------- ! 
     657         DO jk = 1, jpkm1 
     658            e3t_a(:,:,jk) = e3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                      
     659         END DO 
     660      ENDIF 
     661 
     662      IF( ln_vvl_dbg.AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN   ! - ML - test: control prints for debuging 
     663         ! 
     664         zht(:,:) = 0.0_wp 
     665         DO jk = 1, jpkm1 
     666            zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
     667         END DO 
     668         IF( lwp ) WRITE(numout, *) 'kt = ', kt 
     669         IF( lwp ) WRITE(numout, *) 'ncall = ', ncall 
     670         IF( lwp ) WRITE(numout, *) 'll_do_bclinic', ll_do_bclinic  
     671         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     672            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 ) 
     673            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain 
     674            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 
     675         END IF 
     676         ! 
     677         z_tmin = MINVAL( e3t_n(:,:,:)  ) 
    433678         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    434          ! - ML - test: for the moment, stop simulation for too large e3_t variations 
    435          IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 
     679         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3t_n) =', z_tmin 
     680         IF ( z_tmin .LE. 0._wp ) THEN 
    436681            IF( lk_mpp ) THEN 
    437                CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 
    438                CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     682               CALL mpp_minloc(e3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
    439683            ELSE 
    440                ijk_max = MAXLOC( ze3t(:,:,:) ) 
    441                ijk_max(1) = ijk_max(1) + nimpp - 1 
    442                ijk_max(2) = ijk_max(2) + njmpp - 1 
    443                ijk_min = MINLOC( ze3t(:,:,:) ) 
     684               ijk_min = MINLOC( e3t_n(:,:,:)  ) 
    444685               ijk_min(1) = ijk_min(1) + nimpp - 1 
    445686               ijk_min(2) = ijk_min(2) + njmpp - 1 
    446687            ENDIF 
    447688            IF (lwp) THEN 
    448                WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 
    449                WRITE(numout, *) 'at i, j, k=', ijk_max 
    450                WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 
     689               WRITE(numout, *) 'Negative scale factor, e3t_n =', z_tmin 
    451690               WRITE(numout, *) 'at i, j, k=', ijk_min             
    452                CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 
     691               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
     692            ENDIF 
     693         ENDIF          
     694         ! 
     695         z_tmin = MINVAL( e3u_n(:,:,:)) 
     696         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
     697         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3u_n) =', z_tmin 
     698         IF ( z_tmin .LE. 0._wp ) THEN 
     699            IF( lk_mpp ) THEN 
     700               CALL mpp_minloc(e3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     701            ELSE 
     702               ijk_min = MINLOC( e3u_n(:,:,:)  ) 
     703               ijk_min(1) = ijk_min(1) + nimpp - 1 
     704               ijk_min(2) = ijk_min(2) + njmpp - 1 
     705            ENDIF 
     706            IF (lwp) THEN 
     707               WRITE(numout, *) 'Negative scale factor, e3u_n =', z_tmin 
     708               WRITE(numout, *) 'at i, j, k=', ijk_min             
     709               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
     710            ENDIF 
     711         ENDIF  
     712         ! 
     713         z_tmin = MINVAL( e3v_n(:,:,:) ) 
     714         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
     715         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3v_n) =', z_tmin 
     716         IF ( z_tmin .LE. 0._wp ) THEN 
     717            IF( lk_mpp ) THEN 
     718               CALL mpp_minloc(e3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     719            ELSE 
     720               ijk_min = MINLOC( e3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0   ) 
     721               ijk_min(1) = ijk_min(1) + nimpp - 1 
     722               ijk_min(2) = ijk_min(2) + njmpp - 1 
     723            ENDIF 
     724            IF (lwp) THEN 
     725               WRITE(numout, *) 'Negative scale factor, e3v_n =', z_tmin 
     726               WRITE(numout, *) 'at i, j, k=', ijk_min             
     727               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
     728            ENDIF 
     729         ENDIF  
     730         ! 
     731         z_tmin = MINVAL( e3f_n(:,:,:)) 
     732         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
     733         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(e3f_n) =', z_tmin 
     734         IF ( z_tmin .LE. 0._wp ) THEN 
     735            IF( lk_mpp ) THEN 
     736               CALL mpp_minloc(e3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     737            ELSE 
     738               ijk_min = MINLOC( e3f_n(:,:,:) ) 
     739               ijk_min(1) = ijk_min(1) + nimpp - 1 
     740               ijk_min(2) = ijk_min(2) + njmpp - 1 
     741            ENDIF 
     742            IF (lwp) THEN 
     743               WRITE(numout, *) 'Negative scale factor, e3f_n =', z_tmin 
     744               WRITE(numout, *) 'at i, j, k=', ijk_min             
     745               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
    453746            ENDIF 
    454747         ENDIF 
    455          ! - ML - end test 
    456          ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 
    457          tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) ) 
    458          tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 
    459  
    460          ! 
    461          ! "tilda" change in the after scale factor 
    462          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    463          DO jk = 1, jpkm1 
    464             dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 
    465          END DO 
    466          ! III - Barotropic repartition of the sea surface height over the baroclinic profile 
    467          ! ================================================================================== 
    468          ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 
    469          ! - ML - baroclinicity error should be better treated in the future 
    470          !        i.e. locally and not spread over the water column. 
    471          !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 
    472          zht(:,:) = 0. 
    473          DO jk = 1, jpkm1 
    474             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    475          END DO 
    476          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
    477          DO jk = 1, jpkm1 
    478             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    479          END DO 
    480  
    481       ENDIF 
    482  
    483       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
    484       !                                           ! ---baroclinic part--------- ! 
    485          DO jk = 1, jpkm1 
    486             e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    487          END DO 
    488       ENDIF 
    489  
    490       IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging 
    491          ! 
    492          IF( lwp ) WRITE(numout, *) 'kt =', kt 
    493          IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    494             z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 
    495             IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain 
    496             IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 
    497          END IF 
    498748         ! 
    499749         zht(:,:) = 0.0_wp 
    … …  
    503753         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 
    504754         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    505          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 
     755         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn  -SUM(e3t_n))) =', z_tmax 
    506756         ! 
    507757         zht(:,:) = 0.0_wp 
    508758         DO jk = 1, jpkm1 
    509             zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 
    510          END DO 
    511          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 
     759            zht(:,:) = zht(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
     760         END DO 
     761         zwu(:,:) = 0._wp 
     762         DO jj = 1, jpjm1 
     763            DO ji = 1, fs_jpim1   ! vector opt. 
     764               zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e1e2u(ji,jj)                             & 
     765                        &          * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj) ) 
     766            END DO 
     767         END DO 
     768         CALL lbc_lnk( zwu(:,:), 'U', 1._wp ) 
     769         z_tmax = MAXVAL( ssumask(:,:) * ssumask(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) ) 
    512770         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    513          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 
     771         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(e3u_n))) =', z_tmax 
    514772         ! 
    515773         zht(:,:) = 0.0_wp 
    516774         DO jk = 1, jpkm1 
    517             zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 
    518          END DO 
    519          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 
     775            zht(:,:) = zht(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
     776         END DO 
     777         zwv(:,:) = 0._wp 
     778         DO jj = 1, jpjm1 
     779            DO ji = 1, fs_jpim1   ! vector opt. 
     780               zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e1e2v(ji,jj)                             & 
     781                        &          * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1) ) 
     782            END DO 
     783         END DO 
     784         CALL lbc_lnk( zwv(:,:), 'V', 1._wp ) 
     785         z_tmax = MAXVAL( ssvmask(:,:) * ssvmask(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) ) 
    520786         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    521          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 
    522          ! 
    523          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) ) 
     787         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(e3v_n))) =', z_tmax 
     788         ! 
     789         zht(:,:) = 0.0_wp 
     790         DO jk = 1, jpkm1 
     791            DO jj = 1, jpjm1 
     792               zht(:,jj) = zht(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk) 
     793            END DO 
     794         END DO 
     795         zwu(:,:) = 0._wp 
     796         DO jj = 1, jpjm1 
     797            DO ji = 1, fs_jpim1   ! vector opt. 
     798               zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e1e2f(ji,jj)  & 
     799                        &          * (  e1e2t(ji  ,jj) * sshn(ji  ,jj) + e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1) & 
     800                        &             + e1e2t(ji+1,jj) * sshn(ji+1,jj) + e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) ) 
     801            END DO 
     802         END DO 
     803         CALL lbc_lnk( zht(:,:), 'F', 1._wp ) 
     804         CALL lbc_lnk( zwu(:,:), 'F', 1._wp ) 
     805         z_tmax = MAXVAL( fmask(:,:,1) * fmask(:,:,1) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) ) 
    524806         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    525          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 
    526          ! 
    527          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) ) 
    528          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    529          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 
    530          ! 
    531          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) ) 
    532          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    533          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 
     807         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(e3f_n))) =', z_tmax 
     808         ! 
    534809      END IF 
    535810 
    … …  
    699974      !!---------------------------------------------------------------------- 
    700975      ! 
    701       nmet = 0  ! Original method (Surely wrong) 
    702 !     nmet = 1  ! Interface interpolation 
    703 !     nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment 
     976!     nmet = 0  ! Original method (Surely wrong) 
     977     nmet = 2  ! Interface interpolation 
     978!      nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment 
    704979!     Note that we kept surface weighted interpolation for barotropic increment to be compliant 
    705980!     with what is done in surface pressure module. 
    … …  
    711986      END IF 
    712987      ! 
    713       ztap   = 0.0_wp   ! Minimum fraction of T-point thickness at cell interfaces 
    714       zsmall = 1.e-3_wp ! Minimum thickness at U or V points (m) 
     988      ztap   = 0._wp   ! Minimum fraction of T-point thickness at cell interfaces 
     989      zsmall = 1.e-8_wp ! Minimum thickness at U or V points (m) 
    715990      ! 
    716991      IF ( (nmet==1).OR.(nmet==2) ) THEN 
    … …  
    7441019               END DO  
    7451020            ENDIF  
     1021            zw(:,:,:) = (zw(:,:,:) - gdepw_0(:,:,:))*tmask(:,:,:) 
    7461022            ! 
    7471023         END SELECT 
    … …  
    7681044                  ! Correction at last level: 
    7691045                  jkbot = mbku(ji,jj) 
    770                   zdo   = hu_0(ji,jj) 
     1046                  zdo = 0._wp 
    7711047                  DO jk=jkbot,1,-1 
    772                      zup = 0.5_wp * ( zw(ji  ,jj,jk) + zw(ji+1,jj,jk) )  
     1048                     zup = 0.5_wp * ( e1e2t(ji  ,jj)*zw(ji  ,jj,jk) & 
     1049                         &          + e1e2t(ji+1,jj)*zw(ji+1,jj,jk) ) * r1_e1e2u(ji,jj) 
    7731050                     ! 
    7741051                     ! If there is a step, taper bottom interface: 
    775                      IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(zup>zdo)) THEN 
    776                         IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN 
    777                            zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk)) 
    778                         ELSE 
    779                            zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk)) 
    780                         ENDIF 
    781                         zup = MIN(zup, zdo-zmin) 
    782                      ENDIF 
    783                      zup = MIN(zup, zdo-zsmall) 
    784                      pe3_out(ji,jj,jk) = (zdo - zup - e3u_0(ji,jj,jk)) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) 
     1052!                     IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(jk==jkbot)) THEN 
     1053!                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN 
     1054!                           zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk)) 
     1055!                        ELSE 
     1056!                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk)) 
     1057!                        ENDIF 
     1058!                        zup = MIN(zup, zdo-zmin) 
     1059!                     ENDIF 
     1060                     zup = MIN(zup, zdo+e3u_0(ji,jj,jk)-zsmall) 
     1061                     pe3_out(ji,jj,jk) = (zdo - zup) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) 
    7851062                     zdo = zup 
    7861063                  END DO 
    … …  
    7951072                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                   & 
    7961073                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )               &  
    797                            &               / ( hu_0(ji,jj) + 1._wp - umask(ji,jj,1) )              & 
     1074                           &               / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )              & 
    7981075                           &               * 0.5_wp * r1_e1e2u(ji,jj)                              & 
    7991076                           &               * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )        & 
    … …  
    8241101                  ! Correction at last level: 
    8251102                  jkbot = mbkv(ji,jj) 
    826                   zdo   = hv_0(ji,jj) 
     1103                  zdo = 0._wp 
    8271104                  DO jk=jkbot,1,-1 
    828                      zup = 0.5_wp * ( zw(ji,jj  ,jk) + zw(ji,jj+1,jk) )  
     1105                     zup = 0.5_wp * ( e1e2t(ji,jj  ) * zw(ji,jj  ,jk) &  
     1106                         &          + e1e2t(ji,jj+1) * zw(ji,jj+1,jk) ) * r1_e1e2v(ji,jj) 
    8291107                     ! 
    8301108                     ! If there is a step, taper bottom interface: 
    831                      IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN 
    832                         IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN 
    833                            zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk)) 
    834                         ELSE 
    835                            zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk)) 
    836                         ENDIF 
    837                         zup = MIN(zup, zdo-zmin) 
    838                      ENDIF 
    839                      zup = MIN(zup, zdo-zsmall) 
    840                      pe3_out(ji,jj,jk) = (zdo - zup - e3v_0(ji,jj,jk)) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) 
     1109!                     IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN 
     1110!                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN 
     1111!                           zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk)) 
     1112!                        ELSE 
     1113!                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk)) 
     1114!                        ENDIF 
     1115!                        zup = MIN(zup, zdo-zmin) 
     1116!                     ENDIF 
     1117                     zup = MIN(zup, zdo + e3v_0(ji,jj,jk) - zsmall) 
     1118                     pe3_out(ji,jj,jk) = (zdo - zup) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) 
    8411119                     zdo = zup 
    8421120                  END DO 
    … …  
    8511129                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                                       & 
    8521130                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )                                   &  
    853                            &               / ( hv_0(ji,jj) + 1._wp - vmask(ji,jj,1) )                                  & 
     1131                           &               / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )                                  & 
    8541132                           &               * 0.5_wp * r1_e1e2v(ji,jj)                                                  & 
    8551133                                           * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )                            & 
    … …  
    8831161                  ! bottom correction: 
    8841162                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1))  
    885                   zdo = hf_0(ji,jj) 
     1163                  zdo = 0._wp 
    8861164                  DO jk=jkbot,1,-1 
    887                      zup =  0.25_wp * (  zw(ji  ,jj  ,jk)  &  
    888                            &           + zw(ji+1,jj  ,jk)  &  
    889                            &           + zw(ji  ,jj+1,jk)  &  
    890                            &           + zw(ji+1,jj+1,jk)  ) 
     1165                     zup =  0.25_wp * (  e1e2t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  &  
     1166                           &           + e1e2t(ji+1,jj  ) * zw(ji+1,jj  ,jk)  &  
     1167                           &           + e1e2t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  &  
     1168                           &           + e1e2t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  ) *    r1_e1e2f(ji,jj) 
    8911169                     ! 
    8921170                     ! If there is a step, taper bottom interface: 
    893                      IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN 
    894                         IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN 
    895                            IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN 
    896                               zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk)) 
    897                            ELSE 
    898                               zmin = ztap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk)) 
    899                            ENDIF 
    900                         ELSE 
    901                            IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN 
    902                               zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk)) 
    903                            ELSE 
    904                               zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk)) 
    905                            ENDIF 
    906                         ENDIF 
    907                         zup = MIN(zup, zdo-zmin) 
    908                      ENDIF 
    909                      zup = MIN(zup, zdo-zsmall) 
     1171!                     IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN 
     1172!                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN 
     1173!                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN 
     1174!                              zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk)) 
     1175!                           ELSE 
     1176!                              zmin = ztap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk)) 
     1177!                           ENDIF 
     1178!                        ELSE 
     1179!                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN 
     1180!                              zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk)) 
     1181!                           ELSE 
     1182!                              zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk)) 
     1183!                           ENDIF 
     1184!                        ENDIF 
     1185!                        zup = MIN(zup, zdo-zmin) 
     1186!                     ENDIF 
     1187                     zup = MIN(zup, zdo+e3f_0(ji,jj,jk)-zsmall) 
    9101188                     ! 
    911                      pe3_out(ji,jj,jk) = ( zdo - zup - e3f_0(ji,jj,jk) ) &  
     1189                     pe3_out(ji,jj,jk) = ( zdo - zup ) &  
    9121190                                      & *( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) 
    9131191                     zdo = zup 
    … …  
    10061284            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
    10071285            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    1008             id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     1286            id5 = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. ) 
    10091287            !                             ! --------- ! 
    10101288            !                             ! all cases ! 
    … …  
    10681346                  !                       ! ------------ ! 
    10691347                  IF( id5 > 0 ) THEN  ! required array exists 
    1070                      CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios ) 
     1348                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lrxios ) 
    10711349                  ELSE                ! array is missing 
    1072                      hdiv_lf(:,:,:) = 0.0_wp 
     1350                     hdivn_lf(:,:,:) = 0.0_wp 
    10731351                  ENDIF 
    10741352               ENDIF 
    … …  
    11401418               tilde_e3t_b(:,:,:) = 0._wp 
    11411419               tilde_e3t_n(:,:,:) = 0._wp 
    1142                IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 
     1420               IF( ln_vvl_ztilde ) THEN  
     1421                  hdivn_lf(:,:,:) = 0._wp  
     1422                  un_lf(:,:,:) = 0._wp  
     1423                  vn_lf(:,:,:) = 0._wp  
     1424               ENDIF 
    11431425            END IF 
    11441426         ENDIF 
    … …  
    11621444         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    11631445            !                                        ! ------------ ! 
    1164             CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios) 
     1446            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lwxios) 
    11651447         ENDIF 
    11661448         ! 
    … …  
    11791461      !!---------------------------------------------------------------------- 
    11801462      INTEGER ::   ioptio, ios 
    1181       !! 
    1182       NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 
    1183          &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , & 
    1184          &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe 
    1185       !!----------------------------------------------------------------------  
     1463 
     1464      NAMELIST/nam_vvl/ ln_vvl_zstar               , ln_vvl_ztilde                       , & 
     1465                      & ln_vvl_layer               , ln_vvl_ztilde_as_zstar              , & 
     1466                      & ln_vvl_zstar_at_eqtor      , ln_vvl_zstar_on_shelf               , & 
     1467                      & ln_vvl_adv_cn2             , ln_vvl_adv_fct                      , & 
     1468                      & ln_vvl_lap                 , ln_vvl_blp                          , & 
     1469                      & rn_ahe3_lap                , rn_ahe3_blp                         , & 
     1470                      & rn_rst_e3t                 , rn_lf_cutoff                        , & 
     1471                      & ln_vvl_regrid                                                    , & 
     1472                      & ln_vvl_ramp                , rn_day_ramp                         , & 
     1473                      & ln_vvl_dbg   ! not yet implemented: ln_vvl_kepe 
     1474      !!----------------------------------------------------------------------   
    11861475      ! 
    11871476      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :  
    … …  
    11971486         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate' 
    11981487         WRITE(numout,*) '~~~~~~~~~~~' 
    1199          WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate' 
    1200          WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar 
    1201          WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde 
    1202          WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer 
    1203          WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar 
     1488         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate' 
     1489         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar 
     1490         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde 
     1491         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer 
     1492         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar 
     1493         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation' 
     1494         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe 
    12041495         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor 
    1205          WRITE(numout,*) '      !' 
    1206          WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3 
    1207          WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max 
     1496         WRITE(numout,*) '      ztilde on shelves          ln_vvl_zstar_on_shelf  = ', ln_vvl_zstar_on_shelf 
     1497         WRITE(numout,*) '           Namelist nam_vvl : thickness advection scheme' 
     1498         WRITE(numout,*) '              2nd order                  ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2 
     1499         WRITE(numout,*) '              2nd order FCT              ln_vvl_adv_fct = ', ln_vvl_adv_fct                     
     1500         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion scheme' 
     1501         WRITE(numout,*) '              Laplacian                  ln_vvl_lap     = ', ln_vvl_lap 
     1502         WRITE(numout,*) '              Bilaplacian                ln_vvl_blp     = ', ln_vvl_blp 
     1503         WRITE(numout,*) '              Laplacian   coefficient    rn_ahe3_lap    = ', rn_ahe3_lap 
     1504         WRITE(numout,*) '              Bilaplacian coefficient    rn_ahe3_blp    = ', rn_ahe3_blp 
     1505         WRITE(numout,*) '           Namelist nam_vvl : layers regriding' 
     1506         WRITE(numout,*) '              ln_vvl_regrid                             = ', ln_vvl_regrid 
     1507         WRITE(numout,*) '           Namelist nam_vvl : linear ramp at startup' 
     1508         WRITE(numout,*) '              ln_vvl_ramp                               = ', ln_vvl_ramp 
     1509         WRITE(numout,*) '              rn_day_ramp                               = ', rn_day_ramp 
    12081510         IF( ln_vvl_ztilde_as_zstar ) THEN 
    1209             WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) ' 
    1210             WRITE(numout,*) '         ignoring namelist timescale parameters and using:' 
    1211             WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)' 
    1212             WRITE(numout,*) '                         rn_rst_e3t     = 0.e0' 
    1213             WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 
    1214             WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt' 
     1511            WRITE(numout,*) '           ztilde running in zstar emulation mode; ' 
     1512            WRITE(numout,*) '           ignoring namelist timescale parameters and using:' 
     1513            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)' 
     1514            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0' 
     1515            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 
     1516            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt' 
    12151517         ELSE 
    1216             WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t 
    1217             WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff 
     1518            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)' 
     1519            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t 
     1520            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)' 
     1521            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff 
    12181522         ENDIF 
    1219          WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg 
     1523         WRITE(numout,*) '           Namelist nam_vvl : debug prints' 
     1524         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg 
     1525      ENDIF 
     1526      ! 
     1527      IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN  
     1528         ioptio = 0                      ! Choose one advection scheme at most 
     1529         IF( ln_vvl_adv_cn2         )        ioptio = ioptio + 1 
     1530         IF( ln_vvl_adv_fct         )        ioptio = ioptio + 1 
     1531         IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' ) 
    12201532      ENDIF 
    12211533      ! 
    … …  
    12371549      ENDIF 
    12381550      ! 
     1551      ! Use of "shelf horizon depths" should be allowed with s-z coordinates, but we restrict it to zco and zps  
     1552      ! for the time being 
     1553      IF ( ln_sco ) THEN  
     1554        ll_shorizd=.FALSE. 
     1555      ELSE 
     1556        ll_shorizd=.TRUE. 
     1557      ENDIF 
     1558      ! 
    12391559#if defined key_agrif 
    12401560      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' ) 
    … …  
    12431563   END SUBROUTINE dom_vvl_ctl 
    12441564 
     1565   SUBROUTINE dom_vvl_regrid( kt ) 
     1566      !!---------------------------------------------------------------------- 
     1567      !!                ***  ROUTINE dom_vvl_regrid  *** 
     1568      !!                    
     1569      !! ** Purpose :  Ensure "well-behaved" vertical grid 
     1570      !! 
     1571      !! ** Method  :  More or less adapted from references below. 
     1572      !!regrid 
     1573      !! ** Action  :  Ensure that thickness are above a given value, spaced enough 
     1574      !!               and revert to Eulerian coordinates near the bottom.       
     1575      !! 
     1576      !! References :  Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction 
     1577      !!               with a Model Combining Terrain-following and Isentropic 
     1578      !!               coordinates. Part I: Model Description. Monthly Weather Rev., 
     1579      !!               121, 1770-1785. 
     1580      !!               Toy, M., 2011: Incorporating Condensational Heating into a 
     1581      !!               Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic- 
     1582      !!               Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954. 
     1583      !!---------------------------------------------------------------------- 
     1584      !! * Arguments 
     1585      INTEGER, INTENT( in )               :: kt         ! time step 
     1586 
     1587      !! * Local declarations 
     1588      INTEGER                             :: ji, jj, jk ! dummy loop indices 
     1589      LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond 
     1590      LOGICAL                             :: ll_zdiff_cond, ll_blpdiff_cond 
     1591      INTEGER                             :: jkbot 
     1592      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd 
     1593      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj, dzmin_int, dzmin_surf 
     1594      REAL(wp)                            :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn 
     1595      REAL(wp)                            :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot 
     1596      REAL(wp)                            :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim 
     1597      REAL(wp), DIMENSION(jpi,jpj)        :: zdw, zwu, zwv 
     1598      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: zwdw, zwdw_b 
     1599      !!---------------------------------------------------------------------- 
     1600 
     1601      IF( ln_timing )  CALL timing_start('dom_vvl_regrid') 
     1602      ! 
     1603      ! 
     1604      ! Some user defined parameters below: 
     1605      ll_chk_bot2top   = .TRUE. 
     1606      ll_chk_top2bot   = .TRUE. 
     1607      dzmin_int  = 0.1_wp   ! Absolute minimum depth in the interior (in meters) 
     1608      dzmin_surf = 1.0_wp   ! Absolute minimum depth at the surface (in meters) 
     1609      zfrch_stp  = 0.5_wp   ! Maximum fractionnal thickness change in one time step (<= 1.) 
     1610      zfrch_rel  = 0.4_wp   ! Maximum relative thickness change in the vertical (<= 1.) 
     1611      zfrac_bot  = 0.05_wp  ! Fraction of bottom level allowed to change 
     1612      zscal_bot  = 2.0_wp   ! Depth lengthscale 
     1613      ll_zdiff_cond    = .TRUE.  ! Conditionnal vertical diffusion of interfaces 
     1614         zvdiff  = 0.2_wp   ! m 
     1615         zvlim   = 0.5_wp   ! max d2h/dh 
     1616      ll_lapdiff_cond  = .TRUE.  ! Conditionnal Laplacian diffusion of interfaces 
     1617         zhdiff  = 0.01_wp   ! ad. 
     1618         zhlim   = 0.03_wp  ! ad. max lap(z)*e1 
     1619      ll_blpdiff_cond  = .FALSE.  ! Conditionnal Bilaplacian diffusion of interfaces 
     1620         zhdiff2 = 0.2_wp   ! ad. 
     1621!         zhlim2  = 3.e-11_wp  !  max bilap(z) 
     1622         zhlim2  = 0.01_wp  ! ad. max bilap(z)*e1**3 
     1623      ! ---------------------------------------------------------------------------------------  
     1624      ! 
     1625      ! Set arrays determining maximum vertical displacement at the bottom: 
     1626      !-------------------------------------------------------------------- 
     1627      IF ( kt==nit000 ) THEN 
     1628         DO jj = 2, jpjm1 
     1629            DO ji = 2, jpim1 
     1630               jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1)) 
     1631               jk = MIN(jk,mbkt(ji-1,jj-1), mbkt(ji-1,jj+1), mbkt(ji+1,jj+1), mbkt(ji+1,jj-1)) 
     1632               i_int_bot(ji,jj) = jk 
     1633            END DO 
     1634         END DO 
     1635         dsm(:,:) = REAL( i_int_bot(:,:), wp )   ;   CALL lbc_lnk(dsm(:,:),'T',1.) 
     1636         i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 ) 
     1637 
     1638         DO jj = 2, jpjm1 
     1639            DO ji = 2, jpim1 
     1640               zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji  ,jj,1), & 
     1641                          &     ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), & 
     1642                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj  ,1), & 
     1643                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1)  ) 
     1644                zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall ) 
     1645            END DO 
     1646         END DO 
     1647         CALL lbc_lnk( zdw(:,:), 'T', 1. ) 
     1648 
     1649         DO jj = 2, jpjm1 
     1650            DO ji = 2, jpim1 
     1651               dsm(ji,jj) = 1._wp/16._wp * (            zdw(ji-1,jj-1) + zdw(ji+1,jj-1)      & 
     1652                  &                           +         zdw(ji-1,jj+1) + zdw(ji+1,jj+1)      & 
     1653                  &                           + 2._wp*( zdw(ji  ,jj-1) + zdw(ji-1,jj  )      & 
     1654                  &                           +         zdw(ji+1,jj  ) + zdw(ji  ,jj+1) )    & 
     1655                  &                           + 4._wp*  zdw(ji  ,jj  )                       ) 
     1656            END DO 
     1657         END DO          
     1658 
     1659         CALL lbc_lnk( dsm(:,:), 'T', 1. )    
     1660      
     1661         IF (ln_zps) THEN 
     1662            DO jj = 1, jpj 
     1663               DO ji = 1, jpi 
     1664                  jk = i_int_bot(ji,jj) 
     1665                  hsm(ji,jj) = zfrac_bot * e3w_1d(jk) 
     1666!                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj)) 
     1667               END DO 
     1668            END DO 
     1669         ELSE 
     1670            DO jj = 1, jpj 
     1671               DO ji = 1, jpi 
     1672                  jk = i_int_bot(ji,jj) 
     1673                  hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk) 
     1674!                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj)) 
     1675               END DO 
     1676            END DO 
     1677         ENDIF 
     1678      END IF 
     1679 
     1680      ! Provisionnal interface depths: 
     1681      !------------------------------- 
     1682      zwdw(:,:,1) = 0.e0 
     1683      DO jj = 1, jpj 
     1684         DO ji = 1, jpi 
     1685            DO jk = 2, jpk 
     1686               zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + &  
     1687                              & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     1688            END DO 
     1689         END DO 
     1690      END DO 
     1691      ! 
     1692      ! Conditionnal horizontal Laplacian diffusion: 
     1693      !--------------------------------------------- 
     1694      IF ( ll_lapdiff_cond ) THEN 
     1695         ! 
     1696         zwdw_b(:,:,1) = 0._wp 
     1697         DO jj = 1, jpj 
     1698            DO ji = 1, jpi 
     1699               DO jk=2,jpk 
     1700                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + &  
     1701                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     1702               END DO 
     1703            END DO 
     1704         END DO 
     1705         ! 
     1706         DO jk = 2, jpkm1 
     1707            zwu(:,:) = 0._wp 
     1708            zwv(:,:) = 0._wp 
     1709            DO jj = 1, jpjm1 
     1710               DO ji = 1, fs_jpim1   ! vector opt.                   
     1711                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) & 
     1712                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) ) 
     1713                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) &  
     1714                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) ) 
     1715               END DO 
     1716            END DO 
     1717            DO jj = 2, jpjm1 
     1718               DO ji = fs_2, fs_jpim1   ! vector opt. 
     1719                  ztmp1 = ( zwu(ji-1,jj  ) - zwu(ji,jj) & 
     1720                     &  +   zwv(ji  ,jj-1) - zwv(ji,jj) ) * r1_e1e2t(ji,jj) 
     1721                  zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e1e2t(ji,jj)), 0._wp) 
     1722                  ztmp = SIGN(zh2, ztmp1) 
     1723                  zeu2 = zhdiff * e1e2t(ji,jj)*e1e2t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj)) 
     1724                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk) 
     1725               END DO 
     1726            END DO 
     1727         END DO          
     1728         ! 
     1729      ENDIF 
     1730 
     1731      ! Conditionnal horizontal Bilaplacian diffusion: 
     1732      !----------------------------------------------- 
     1733      IF ( ll_blpdiff_cond ) THEN 
     1734         ! 
     1735         zwdw_b(:,:,1) = 0._wp 
     1736         DO jj = 1, jpj 
     1737            DO ji = 1, jpi 
     1738               DO jk = 2,jpkm1 
     1739                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + &  
     1740                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     1741               END DO 
     1742            END DO 
     1743         END DO 
     1744         ! 
     1745         DO jk = 2, jpkm1 
     1746            zwu(:,:) = 0._wp 
     1747            zwv(:,:) = 0._wp 
     1748            DO jj = 1, jpjm1 
     1749               DO ji = 1, fs_jpim1   ! vector opt.                   
     1750                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) & 
     1751                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) ) 
     1752                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) &  
     1753                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) ) 
     1754               END DO 
     1755            END DO 
     1756            DO jj = 2, jpjm1 
     1757               DO ji = fs_2, fs_jpim1   ! vector opt. 
     1758                  zwdw_b(ji,jj,jk) = -( (zwu(ji-1,jj  ) - zwu(ji,jj)) & 
     1759                                &  +    (zwv(ji  ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj) 
     1760               END DO 
     1761            END DO 
     1762         END DO    
     1763         ! 
     1764         CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )       
     1765         ! 
     1766         DO jk = 2, jpkm1 
     1767            zwu(:,:) = 0._wp 
     1768            zwv(:,:) = 0._wp 
     1769            DO jj = 1, jpjm1 
     1770               DO ji = 1, fs_jpim1   ! vector opt.                   
     1771                  zwu(ji,jj) =  umask(ji,jj,jk) * e2_e1u(ji,jj) & 
     1772                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) ) 
     1773                  zwv(ji,jj) =  vmask(ji,jj,jk) * e1_e2v(ji,jj) &  
     1774                             &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) ) 
     1775               END DO 
     1776            END DO 
     1777            DO jj = 2, jpjm1 
     1778               DO ji = fs_2, fs_jpim1   ! vector opt. 
     1779                  ztmp1 = ( (zwu(ji-1,jj  ) - zwv(ji,jj)) & 
     1780                     &  +   (zwv(ji  ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj) 
     1781                  zh2 = MAX(abs(ztmp1)-zhlim2*SQRT(r1_e1e2t(ji,jj))*r1_e1e2t(ji,jj), 0._wp) 
     1782                  ztmp = SIGN(zh2, ztmp1) 
     1783                  zeu2 = zhdiff2 * e1e2t(ji,jj)*e1e2t(ji,jj) / 16._wp 
     1784                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk) 
     1785               END DO 
     1786            END DO 
     1787         END DO    
     1788         ! 
     1789      ENDIF 
     1790 
     1791      ! Conditionnal vertical diffusion: 
     1792      !--------------------------------- 
     1793      IF ( ll_zdiff_cond ) THEN 
     1794         DO jk = 2, jpkm1 
     1795            DO jj = 2, jpjm1 
     1796               DO ji = fs_2, fs_jpim1   ! vector opt.     
     1797                  ztmp  = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) &  
     1798                            -(tilde_e3t_b(ji,jj,jk  )+e3t_0(ji,jj,jk  ))*tmask(ji,jj,jk  ) )  
     1799                  ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) &  
     1800                        &           +tilde_e3t_b(ji,jj,jk  ) + e3t_0(ji,jj,jk  ) ) 
     1801                  zh2   = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp) 
     1802                  ztmp  = SIGN(zh2, ztmp) 
     1803                  IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0 
     1804                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk) 
     1805               END DO 
     1806            END DO 
     1807         END DO 
     1808      ENDIF 
     1809      ! 
     1810      ! Check grid from the bottom to the surface 
     1811      !------------------------------------------ 
     1812      IF ( ll_chk_bot2top ) THEN 
     1813         DO jj = 2, jpjm1 
     1814            DO ji = 2, jpim1 
     1815               jkbot = mbkt(ji,jj)                    
     1816               DO jk = jkbot,2,-1 
     1817                  ! 
     1818                  zh_0   = e3t_0(ji,jj,jk) 
     1819                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)) 
     1820                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
     1821                  zh_min = MIN(zh_0/3._wp, dzmin_int) 
     1822                  !  
     1823                  ! Set maximum and minimum vertical excursions    
     1824                  ztmph = hsm(ji,jj) 
     1825                  ztmpd = dsm(ji,jj) 
     1826                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)))/ztmpd) 
     1827!                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd) 
     1828                  zh2 = MAX(zh2,0.001_wp) ! Extend tolerance a bit for stability reasons (to be explored) 
     1829                  zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 ) 
     1830                  zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff) 
     1831                  zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 ) 
     1832                  zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff) 
     1833                  ! 
     1834                  ! New layer thickness: 
     1835                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
     1836                  ! 
     1837                  ! Ensure minimum layer thickness:                   
     1838!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new) 
     1839                  zh_new = cush(zh_new, zh_min) 
     1840                  ! 
     1841                  ! Final flux: 
     1842                  zdiff = (zh_new - zh_old) * tmask(ji,jj,jk) 
     1843                  ! 
     1844                  ! Limit thickness change in 1 time step:  
     1845                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 
     1846                  zdiff = SIGN(ztmp, zh_new - zh_old) 
     1847                  zh_new = zdiff + zh_old 
     1848                  ! 
     1849                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk+1) - zh_new        
     1850               END DO     
     1851            END DO  
     1852         END DO 
     1853      END IF     
     1854      ! 
     1855      ! Check grid from the surface to the bottom  
     1856      !------------------------------------------  
     1857      IF ( ll_chk_top2bot ) THEN       
     1858         DO jj = 2, jpjm1 
     1859            DO ji = 2, jpim1 
     1860               jkbot = mbkt(ji,jj)    
     1861               DO jk = 1, jkbot-1 
     1862                  ! 
     1863                  zh_0   = e3t_0(ji,jj,jk) 
     1864                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)) 
     1865                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
     1866                  zh_min = MIN(zh_0/3._wp, dzmin_int) 
     1867                  ! 
     1868                  zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf) 
     1869                  ! 
     1870                  ! New layer thickness: 
     1871                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
     1872                  ! 
     1873                  ! Ensure minimum layer thickness: 
     1874!                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new) 
     1875                  zh_new = cush(zh_new, zh_min) 
     1876                  ! 
     1877                  ! Final flux: 
     1878                  zdiff = (zh_new -zh_old) * tmask(ji,jj,jk) 
     1879                  !  
     1880                  ! Limit flux:                  
     1881!                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 
     1882!                  zdiff = SIGN(ztmp, zh_new - zh_old) 
     1883                  zh_new = zdiff + zh_old 
     1884                  ! 
     1885                  zwdw(ji,jj,jk+1) = zwdw(ji,jj,jk) + zh_new 
     1886               END DO 
     1887               !                
     1888            END DO 
     1889         END DO 
     1890      ENDIF 
     1891      ! 
     1892      DO jj = 2, jpjm1 
     1893         DO ji = 2, jpim1 
     1894            DO jk = 1, jpkm1 
     1895               tilde_e3t_a(ji,jj,jk) =  (zwdw(ji,jj,jk+1)-zwdw(ji,jj,jk)-e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
     1896            END DO 
     1897         END DO 
     1898      END DO 
     1899      !  
     1900      ! 
     1901      IF( ln_timing )  CALL timing_stop('dom_vvl_regrid') 
     1902      ! 
     1903   END SUBROUTINE dom_vvl_regrid 
     1904 
     1905   FUNCTION cush(hin, hmin)  RESULT(hout) 
     1906      !!---------------------------------------------------------------------- 
     1907      !!                 ***  FUNCTION cush  *** 
     1908      !! 
     1909      !! ** Purpose :    
     1910      !! 
     1911      !! ** Method  : 
     1912      !! 
     1913      !!---------------------------------------------------------------------- 
     1914      IMPLICIT NONE 
     1915      REAL(wp), INTENT(in) ::  hin, hmin 
     1916      REAL(wp)             ::  hout, zx, zh_cri 
     1917      !!---------------------------------------------------------------------- 
     1918      zh_cri = 3._wp * hmin 
     1919      ! 
     1920      IF ( hin<=0._wp ) THEN 
     1921         hout = hmin 
     1922      ! 
     1923      ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN 
     1924         zx = hin/zh_cri 
     1925         hout = hmin * (1._wp + zx + zx*zx)       
     1926      ! 
     1927      ELSEIF ( hin>zh_cri ) THEN 
     1928         hout = hin 
     1929      ! 
     1930      ENDIF 
     1931      ! 
     1932   END FUNCTION cush 
     1933 
     1934   FUNCTION cush_max(hin, hmax)  RESULT(hout) 
     1935      !!---------------------------------------------------------------------- 
     1936      !!                 ***  FUNCTION cush  *** 
     1937      !! 
     1938      !! ** Purpose :    
     1939      !! 
     1940      !! ** Method  : 
     1941      !! 
     1942      !!---------------------------------------------------------------------- 
     1943      IMPLICIT NONE 
     1944      REAL(wp), INTENT(in) ::  hin, hmax 
     1945      REAL(wp)             ::  hout, hmin, zx, zh_cri 
     1946      !!---------------------------------------------------------------------- 
     1947      hmin = 0.1_wp * hmax 
     1948      zh_cri = 3._wp * hmin 
     1949      ! 
     1950      IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN 
     1951         zx = (hmax-hin)/zh_cri 
     1952         hout = hmax - hmin * (1._wp + zx + zx*zx) 
     1953         ! 
     1954      ELSEIF ( hin>(hmax-zh_cri) ) THEN 
     1955         hout = hmax - hmin 
     1956         ! 
     1957      ELSE 
     1958         hout = hin 
     1959         ! 
     1960      ENDIF 
     1961      ! 
     1962   END FUNCTION cush_max 
     1963 
     1964   SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin ) 
     1965      !!---------------------------------------------------------------------- 
     1966      !!                  ***  ROUTINE dom_vvl_adv_fct  *** 
     1967      !!  
     1968      !! **  Purpose :  Do thickness advection 
     1969      !! 
     1970      !! **  Method  :  FCT scheme to ensure positivity  
     1971      !! 
     1972      !! **  Action  : - Update pta thickness tendency and diffusive fluxes 
     1973      !!               - this is the total trend, hence it does include sea level motions 
     1974      !!---------------------------------------------------------------------- 
     1975      ! 
     1976      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     1977      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta    ! thickness baroclinic trend  
     1978      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input velocities 
     1979      ! 
     1980      INTEGER  ::   ji, jj, jk, ib, ib_bdy               ! dummy loop indices   
     1981      INTEGER  ::   ikbu, ikbv, ibot 
     1982      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar 
     1983      REAL(wp) ::   zdi, zdj, zmin           !   -      - 
     1984      REAL(wp) ::   zfp_ui, zfp_vj           !   -      - 
     1985      REAL(wp) ::   zfm_ui, zfm_vj           !   -      - 
     1986      REAL(wp) ::   zfp_hi, zfp_hj           !   -      - 
     1987      REAL(wp) ::   zfm_hi, zfm_hj           !   -      - 
     1988      REAL(wp) ::   ztout,  ztin, zfac       !   -      - 
     1989      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwi 
     1990      !!---------------------------------------------------------------------- 
     1991      ! 
     1992      IF( ln_timing )  CALL timing_start('dom_vvl_adv_fct') 
     1993      ! 
     1994      ! 
     1995      ! 1. Initializations 
     1996      ! ------------------ 
     1997      ! 
     1998      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     1999         z2dtt =  rdt 
     2000      ELSE 
     2001         z2dtt = 2.0_wp * rdt 
     2002      ENDIF 
     2003      ! 
     2004      zwi(:,:,:) = 0.e0 
     2005      zwx(:,:,:) = 0.e0  
     2006      zwy(:,:,:) = 0.e0 
     2007      ! 
     2008      ! 
     2009      ! 2. upstream advection with initial mass fluxes & intermediate update 
     2010      ! -------------------------------------------------------------------- 
     2011      IF ( ll_shorizd ) THEN 
     2012         DO jk = 1, jpkm1 
     2013            DO jj = 1, jpjm1 
     2014               DO ji = 1, fs_jpim1   ! vector opt. 
     2015                  ! 
     2016                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp) 
     2017                  zfp_hi = MIN(zfp_hi, e3t_b(ji  ,jj  ,jk)) 
     2018                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )  
     2019                  ! 
     2020                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp) 
     2021                  zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj  ,jk)) 
     2022                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )  
     2023                  !  
     2024                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp) 
     2025                  zfp_hj = MIN(zfp_hj, e3t_b(ji  ,jj  ,jk)) 
     2026                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )  
     2027                  ! 
     2028                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp) 
     2029                  zfm_hj = MIN(zfm_hj, e3t_b(ji  ,jj+1,jk)) 
     2030                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 
     2031                  ! 
     2032                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 
     2033                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 
     2034                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 
     2035                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 
     2036                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2037                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2038               END DO 
     2039            END DO 
     2040         END DO 
     2041      ELSE 
     2042         DO jk = 1, jpkm1 
     2043            DO jj = 1, jpjm1 
     2044               DO ji = 1, fs_jpim1   ! vector opt. 
     2045                  ! 
     2046                  zfp_hi = e3t_b(ji  ,jj  ,jk) 
     2047                  zfm_hi = e3t_b(ji+1,jj  ,jk)              
     2048                  zfp_hj = e3t_b(ji  ,jj  ,jk) 
     2049                  zfm_hj = e3t_b(ji  ,jj+1,jk) 
     2050                  ! 
     2051                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 
     2052                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 
     2053                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 
     2054                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 
     2055                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2056                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2057               END DO 
     2058            END DO 
     2059         END DO 
     2060      ENDIF 
     2061 
     2062      ! total advective trend 
     2063      DO jk = 1, jpkm1 
     2064         DO jj = 2, jpjm1 
     2065            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2066               zbtr = r1_e1e2t(ji,jj) 
     2067               ! total intermediate advective trends 
     2068               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )  & 
     2069                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )  ) 
     2070               ! 
     2071               ! update and guess with monotonic sheme 
     2072               pta(ji,jj,jk) =   pta(ji,jj,jk) + ztra 
     2073               zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     2074            END DO 
     2075         END DO 
     2076      END DO 
     2077 
     2078      CALL lbc_lnk( zwi, 'T', 1. )   
     2079 
     2080      IF ( ln_bdy ) THEN 
     2081         DO ib_bdy=1, nb_bdy 
     2082            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 
     2083               ji = idx_bdy(ib_bdy)%nbi(ib,1) 
     2084               jj = idx_bdy(ib_bdy)%nbj(ib,1) 
     2085               DO jk = 1, jpkm1   
     2086                  zwi(ji,jj,jk) = e3t_a(ji,jj,jk) 
     2087               END DO 
     2088            END DO  
     2089         END DO 
     2090      ENDIF 
     2091 
     2092      IF ( ln_vvl_dbg ) THEN 
     2093         zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 )  
     2094         IF( lk_mpp )   CALL mpp_min( zmin ) 
     2095         IF( zmin < 0._wp) THEN 
     2096            IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here') 
     2097            IF(lwp) WRITE(numout,*) zmin 
     2098         ENDIF 
     2099      ENDIF 
     2100 
     2101      ! 3. antidiffusive flux : high order minus low order 
     2102      ! -------------------------------------------------- 
     2103      ! antidiffusive flux on i and j 
     2104      DO jk = 1, jpkm1 
     2105         DO jj = 1, jpjm1 
     2106            DO ji = 1, fs_jpim1   ! vector opt. 
     2107               zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * e3u_n(ji,jj,jk) & 
     2108                                & - zwx(ji,jj,jk)) * umask(ji,jj,jk) 
     2109               zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * e3v_n(ji,jj,jk) &  
     2110                                & - zwy(ji,jj,jk)) * vmask(ji,jj,jk) 
     2111               ! 
     2112               ! Update advective fluxes 
     2113               un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk) 
     2114               vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk) 
     2115            END DO 
     2116         END DO 
     2117      END DO 
     2118       
     2119      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions 
     2120 
     2121      ! 4. monotonicity algorithm 
     2122      ! ------------------------- 
     2123      CALL nonosc_2d( e3t_b(:,:,:), zwx, zwy, zwi, z2dtt ) 
     2124 
     2125      ! 5. final trend with corrected fluxes 
     2126      ! ------------------------------------ 
     2127      ! 
     2128      ! Update advective fluxes 
     2129      un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:) 
     2130      vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:) 
     2131      ! 
     2132      DO jk = 1, jpkm1 
     2133         DO jj = 2, jpjm1 
     2134            DO ji = fs_2, fs_jpim1   ! vector opt.   
     2135               ! 
     2136               zbtr = r1_e1e2t(ji,jj) 
     2137               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     2138                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
     2139               ! add them to the general tracer trends 
     2140               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 
     2141            END DO 
     2142         END DO 
     2143      END DO 
     2144      ! 
     2145      IF( ln_timing )  CALL timing_stop('dom_vvl_adv_fct') 
     2146      ! 
     2147   END SUBROUTINE dom_vvl_adv_fct 
     2148 
     2149   SUBROUTINE dom_vvl_ups_cor( kt, pta, uin, vin )  
     2150      !!---------------------------------------------------------------------- 
     2151      !!                  ***  ROUTINE dom_vvl_adv_fct  *** 
     2152      !!  
     2153      !! **  Purpose :  Correct for addionnal barotropic fluxes  
     2154      !!                in the upstream direction 
     2155      !! 
     2156      !! **  Method  :    
     2157      !! 
     2158      !! **  Action  : - Update diffusive fluxes uin, vin 
     2159      !!               - Remove divergence from thickness tendency 
     2160      !!---------------------------------------------------------------------- 
     2161      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     2162      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta       ! thickness baroclinic trend  
     2163      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input fluxes 
     2164      INTEGER  ::   ji, jj, jk               ! dummy loop indices   
     2165      INTEGER  ::   ikbu, ikbv, ibot 
     2166      REAL(wp) ::   zbtr, ztra               ! local scalar 
     2167      REAL(wp) ::   zdi, zdj                 !   -      - 
     2168      REAL(wp) ::   zfp_hi, zfp_hj           !   -      - 
     2169      REAL(wp) ::   zfm_hi, zfm_hj           !   -      - 
     2170      REAL(wp) ::   zfp_ui, zfp_vj           !   -      - 
     2171      REAL(wp) ::   zfm_ui, zfm_vj           !   -      - 
     2172      REAL(wp), DIMENSION(jpi,jpj) :: zbu, zbv, zhu_b, zhv_b 
     2173      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy 
     2174      !!---------------------------------------------------------------------- 
     2175      ! 
     2176      IF( ln_timing )  CALL timing_start('dom_vvl_ups_cor') 
     2177      ! 
     2178      ! Compute barotropic flux difference: 
     2179      zbu(:,:) = 0.e0 
     2180      zbv(:,:) = 0.e0 
     2181      DO jj = 1, jpj 
     2182         DO ji = 1, jpi   ! vector opt. 
     2183            DO jk = 1, jpkm1 
     2184               zbu(ji,jj) = zbu(ji,jj) - uin(ji,jj,jk) * umask(ji,jj,jk) 
     2185               zbv(ji,jj) = zbv(ji,jj) - vin(ji,jj,jk) * vmask(ji,jj,jk) 
     2186            END DO 
     2187         END DO 
     2188      ENDDO 
     2189 
     2190      ! Compute upstream depths: 
     2191      zhu_b(:,:) = 0.e0 
     2192      zhv_b(:,:) = 0.e0 
     2193 
     2194      IF ( ll_shorizd ) THEN 
     2195         ! Correct bottom value 
     2196         ! considering "shelf horizon depth" 
     2197         DO jj = 1, jpjm1 
     2198            DO ji = 1, jpim1   ! vector opt. 
     2199               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 
     2200               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 
     2201               DO jk=1, jpkm1 
     2202                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp) 
     2203                  zfp_hi = MIN(zfp_hi, e3t_b(ji  ,jj  ,jk)) 
     2204                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )  
     2205                  ! 
     2206                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp) 
     2207                  zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj  ,jk)) 
     2208                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 
     2209                  !  
     2210                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp) 
     2211                  zfp_hj = MIN(zfp_hj, e3t_b(ji  ,jj  ,jk)) 
     2212                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )  
     2213                  ! 
     2214                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp) 
     2215                  zfm_hj = MIN(zfm_hj, e3t_b(ji  ,jj+1,jk)) 
     2216                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 
     2217                  ! 
     2218                  zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi & 
     2219                               &                 + (1._wp-zdi) * zfm_hi & 
     2220                               &                ) * umask(ji,jj,jk) 
     2221                  zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj & 
     2222                               &                 + (1._wp-zdj) * zfm_hj & 
     2223                               &                ) * vmask(ji,jj,jk) 
     2224               END DO 
     2225            END DO 
     2226         END DO 
     2227      ELSE 
     2228         DO jj = 1, jpjm1 
     2229            DO ji = 1, jpim1   ! vector opt. 
     2230               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))  
     2231               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 
     2232               DO jk = 1, jpkm1 
     2233                  zfp_hi = e3t_b(ji  ,jj  ,jk) 
     2234                  zfm_hi = e3t_b(ji+1,jj  ,jk) 
     2235                  zfp_hj = e3t_b(ji  ,jj  ,jk) 
     2236                  zfm_hj = e3t_b(ji  ,jj+1,jk) 
     2237                  ! 
     2238                  zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi & 
     2239                               &                  + (1._wp-zdi) * zfm_hi & 
     2240                               &                ) * umask(ji,jj,jk) 
     2241                  ! 
     2242                  zhv_b(ji,jj) = zhv_b(ji,jj) + (          zdj  * zfp_hj & 
     2243                               &                  + (1._wp-zdj) * zfm_hj & 
     2244                               &                ) * vmask(ji,jj,jk) 
     2245               END DO 
     2246            END DO 
     2247         END DO 
     2248      ENDIF 
     2249 
     2250      ! Corrective barotropic velocity (times hor. scale factor) 
     2251      zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1)) 
     2252      zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1)) 
     2253 
     2254      CALL lbc_lnk( zbu(:,:), 'U', -1. ) 
     2255      CALL lbc_lnk( zbv(:,:), 'V', -1. ) 
     2256       
     2257      ! Set corrective fluxes in upstream direction: 
     2258      ! 
     2259      zwx(:,:,:) = 0.e0 
     2260      zwy(:,:,:) = 0.e0 
     2261      IF ( ll_shorizd ) THEN 
     2262         DO jj = 1, jpjm1 
     2263            DO ji = 1, fs_jpim1   ! vector opt. 
     2264               ! upstream scheme 
     2265               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 
     2266               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 
     2267               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 
     2268               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 
     2269               DO jk = 1, jpkm1 
     2270                  zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp) 
     2271                  zfp_hi = MIN(e3t_b(ji  ,jj  ,jk), zfp_hi) 
     2272                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )  
     2273                  ! 
     2274                  zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj  ,jk), 0._wp) 
     2275                  zfm_hi = MIN(e3t_b(ji+1,jj  ,jk), zfm_hi) 
     2276                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )  
     2277                  ! 
     2278                  zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj  ,jk), 0._wp) 
     2279                  zfp_hj = MIN(e3t_b(ji  ,jj  ,jk), zfp_hj)  
     2280                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )  
     2281                  ! 
     2282                  zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji  ,jj+1,jk), 0._wp) 
     2283                  zfm_hj = MIN(e3t_b(ji  ,jj+1,jk), zfm_hj) 
     2284                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )  
     2285                  ! 
     2286                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2287                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2288               END DO 
     2289            END DO 
     2290         END DO 
     2291      ELSE 
     2292         DO jj = 1, jpjm1 
     2293            DO ji = 1, fs_jpim1   ! vector opt. 
     2294               ! upstream scheme 
     2295               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 
     2296               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 
     2297               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 
     2298               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 
     2299               DO jk = 1, jpkm1 
     2300                  zfp_hi = e3t_b(ji  ,jj  ,jk) 
     2301                  zfm_hi = e3t_b(ji+1,jj  ,jk) 
     2302                  zfp_hj = e3t_b(ji  ,jj  ,jk) 
     2303                  zfm_hj = e3t_b(ji  ,jj+1,jk) 
     2304                  ! 
     2305                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2306                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2307               END DO 
     2308            END DO 
     2309         END DO 
     2310      ENDIF 
     2311      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions 
     2312 
     2313      uin(:,:,:) = uin(:,:,:) + zwx(:,:,:) 
     2314      vin(:,:,:) = vin(:,:,:) + zwy(:,:,:) 
     2315      ! 
     2316      ! Update trend with corrective fluxes: 
     2317      DO jk = 1, jpkm1 
     2318         DO jj = 2, jpjm1 
     2319            DO ji = fs_2, fs_jpim1   ! vector opt.   
     2320               ! 
     2321               zbtr = r1_e1e2t(ji,jj) 
     2322               ! total advective trends 
     2323               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     2324                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
     2325               ! add them to the general tracer trends 
     2326               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 
     2327            END DO 
     2328         END DO 
     2329      END DO 
     2330      ! 
     2331      IF( ln_timing )  CALL timing_stop('dom_vvl_ups_cor') 
     2332      ! 
     2333   END SUBROUTINE dom_vvl_ups_cor 
     2334 
     2335   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt ) 
     2336      !!--------------------------------------------------------------------- 
     2337      !!                    ***  ROUTINE nonosc_2d  *** 
     2338      !!      
     2339      !! **  Purpose :   compute monotonic thickness fluxes from the upstream  
     2340      !!       scheme and the before field by a nonoscillatory algorithm  
     2341      !! 
     2342      !! **  Method  :   ... ??? 
     2343      !!       warning : pbef and paft must be masked, but the boundaries 
     2344      !!       conditions on the fluxes are not necessary zalezak (1979) 
     2345      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
     2346      !!       in-space based differencing for fluid 
     2347      !!---------------------------------------------------------------------- 
     2348      ! 
     2349      !!---------------------------------------------------------------------- 
     2350      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
     2351      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field 
     2352      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb        ! monotonic fluxes in the 3 directions 
     2353      ! 
     2354      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     2355      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars 
     2356      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
     2357      REAL(wp) ::   zupip1, zupim1, zupjp1, zupjm1, zupb, zupa 
     2358      REAL(wp) ::   zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa 
     2359      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo 
     2360      !!---------------------------------------------------------------------- 
     2361      ! 
     2362      IF( ln_timing )  CALL timing_start('nonosc2') 
     2363      ! 
     2364      zbig  = 1.e+40_wp 
     2365      zrtrn = 1.e-15_wp 
     2366      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp 
     2367 
     2368 
     2369      ! Search local extrema 
     2370      ! -------------------- 
     2371      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
     2372      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   & 
     2373         &        paft * tmask - zbig * ( 1.e0 - tmask )  ) 
     2374      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   & 
     2375         &        paft * tmask + zbig * ( 1.e0 - tmask )  ) 
     2376 
     2377      DO jk = 1, jpkm1 
     2378         z2dtt = p2dt 
     2379         DO jj = 2, jpjm1 
     2380            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2381 
     2382               ! search maximum in neighbourhood 
     2383               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
     2384                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   & 
     2385                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  )) 
     2386 
     2387               ! search minimum in neighbourhood 
     2388               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   & 
     2389                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   & 
     2390                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  )) 
     2391 
     2392               ! positive part of the flux 
     2393               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   & 
     2394                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )    
     2395 
     2396               ! negative part of the flux 
     2397               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   & 
     2398                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) ) 
     2399 
     2400               ! up & down beta terms 
     2401               zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt 
     2402               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
     2403               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
     2404            END DO 
     2405         END DO 
     2406      END DO 
     2407 
     2408      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
     2409 
     2410      ! 3. monotonic flux in the i & j direction (paa & pbb) 
     2411      ! ---------------------------------------- 
     2412      DO jk = 1, jpkm1 
     2413         DO jj = 2, jpjm1 
     2414            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2415               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
     2416               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
     2417               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) ) 
     2418               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu ) 
     2419 
     2420               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
     2421               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
     2422               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 
     2423               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv ) 
     2424            END DO 
     2425         END DO 
     2426      END DO 
     2427      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign) 
     2428      ! 
     2429      IF( ln_timing )  CALL timing_stop('nonosc2') 
     2430      ! 
     2431   END SUBROUTINE nonosc_2d 
     2432 
    12452433   !!====================================================================== 
    12462434END MODULE domvvl 

  • NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DYN/dynnxt.F90

    r9598 r10116  
    9595      ! 
    9696      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    97       INTEGER  ::   ikt          ! local integers 
     97      INTEGER  ::   ikt, jkbot   ! local integers 
    9898      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars 
    9999      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
     100      REAL(wp) ::   zhcri, zmsku, zwgtu, zmskv, zwgtv  !   -      - 
    100101      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
    101102      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva  
    … …  
    136137            END DO   
    137138         ENDIF 
     139      ENDIF 
     140 
     141      IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN  
     142!         ! Duplicate baroclinic velocities from above in massless layers 
     143!         ! Correction to handle correct barotropic velocities right after 
     144         zhcri = 0.1_wp 
     145         DO ji = 2, jpim1 
     146            DO jj= 2, jpjm1 
     147               jkbot=jpkm1 
     148               DO jk=jpkm1,2,-1 
     149                  IF (e3u_a(ji,jj,jk)>zhcri) jkbot = jk+1 
     150               END DO 
     151               DO jk=jkbot,jpkm1 
     152                  zwgtu = MIN(zhcri, e3u_a(ji,jj,jk)) 
     153                  zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_a(ji,jj,jk) - 0.01_wp) ) 
     154                  ua(ji,jj,jk) = zmsku * (zwgtu * ua(ji,jj,jk) + (zhcri - zwgtu) * ua(ji,jj,jk-1)*zwgtu/(zwgtu+0.01_wp)) / zhcri 
     155               END DO 
     156               jkbot=jpkm1 
     157               DO jk=jpkm1,2,-1 
     158                  IF (e3v_a(ji,jj,jk)>zhcri) jkbot = jk+1 
     159               END DO 
     160               DO jk=jkbot,jpkm1 
     161                  zwgtv = MIN(zhcri, e3v_a(ji,jj,jk)) 
     162                  zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_a(ji,jj,jk) - 0.01_wp) ) 
     163                  va(ji,jj,jk) = zmskv * (zwgtv * va(ji,jj,jk) + (zhcri - zwgtv) * va(ji,jj,jk-1)*zwgtv/(zwgtv+0.01_wp)) / zhcri 
     164               END DO 
     165            END DO 
     166         END DO 
    138167      ENDIF 
    139168 
    … …  
    291320                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    292321                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
     322                        zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, ze3u_f(ji,jj,jk) - 0.01_wp) ) 
     323                        zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, ze3v_f(ji,jj,jk) - 0.01_wp) ) 
    293324                        ! 
    294                         ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity 
    295                         vb(ji,jj,jk) = zvf 
     325                        ub(ji,jj,jk) = zmsku * zuf             ! ub <-- filtered velocity 
     326                        vb(ji,jj,jk) = zmskv * zvf 
    296327                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua 
    297328                        vn(ji,jj,jk) = va(ji,jj,jk) 

  • NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/SBC/sbcwave.F90

    r9966 r10116  
    2121   USE zdf_oce,  ONLY : ln_zdfswm 
    2222   USE bdy_oce        ! open boundary condition variables 
    23    USE domvvl         ! domain: variable volume layers 
    2423   ! 
    2524   USE iom            ! I/O manager library 
Note: See TracChangeset for help on using the changeset viewer.