Changeset 10116


Ignore:
Timestamp:
2018-09-12T17:39:01+02:00 (22 months 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.