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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9353 for branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 – NEMO

Ignore:
Timestamp:
2018-02-23T16:52:00+01:00 (6 years ago)
Author:
jchanut
Message:

Tilde coordinate update: first draft

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r6348 r9353  
    99   !!                                          vvl option includes z_star and z_tilde coordinates 
    1010   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
     11   !!                 !  2018-01  (J. Chanut) improve ztilde robustness 
    1112   !!---------------------------------------------------------------------- 
    1213   !!   'key_vvl'                              variable volume 
     
    3334   USE wrk_nemo        ! Memory allocation 
    3435   USE timing          ! Timing 
     36   USE bdy_oce         ! ocean open boundary conditions 
    3537 
    3638   IMPLICIT NONE 
     
    4850   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.             ! ztilde vertical coordinate 
    4951   LOGICAL , PUBLIC                                      :: ln_vvl_layer = .FALSE.              ! level  vertical coordinate 
    50    LOGICAL , PUBLIC                                      :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate 
    51    LOGICAL , PUBLIC                                      :: ln_vvl_zstar_at_eqtor = .FALSE.     ! ztilde vertical coordinate 
    52    LOGICAL , PUBLIC                                      :: ln_vvl_kepe = .FALSE.               ! kinetic/potential energy transfer 
     52   LOGICAL                                               :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate 
     53   LOGICAL                                               :: ln_vvl_zstar_at_eqtor = .FALSE.     ! revert to zstar at equator 
     54   LOGICAL                                               :: ln_vvl_zstar_on_shelf =.FALSE.      ! revert to zstar on shelves 
     55   LOGICAL                                               :: ln_vvl_kepe =.FALSE.                ! kinetic/potential energy transfer 
     56   LOGICAL                                               :: ln_vvl_adv_fct =.FALSE.             ! Centred thickness advection 
     57   LOGICAL                                               :: ln_vvl_adv_cn2 =.TRUE.              ! FCT thickness advection 
     58   LOGICAL                                               :: ln_vvl_dbg = .FALSE.                ! debug control prints 
     59   LOGICAL                                               :: ln_vvl_lap                          ! Laplacian thickness diffusion 
     60   LOGICAL                                               :: ln_vvl_blp                          ! Bilaplacian thickness diffusion  
     61   LOGICAL                                               :: ln_vvl_regrid                       ! ensure layer separation  
    5362   !                                                                                            ! conservation: not used yet 
    54    REAL(wp)                                              :: rn_ahe3                   ! thickness diffusion coefficient 
     63   REAL(wp)                                              :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian) 
     64   REAL(wp)                                              :: rn_ahe3_blp               ! thickness diffusion coefficient (Bilaplacian) 
    5565   REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days] 
    5666   REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days] 
    5767   REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation 
    58    LOGICAL , PUBLIC                                      :: ln_vvl_dbg = .FALSE.      ! debug control prints 
    5968 
    6069   !! * Module variables 
    61    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                       ! thickness diffusion transport 
    62    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence 
    63    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors 
    64    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors 
    65    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors 
    66    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence 
     70   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport 
     71   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn_lf,  un_lf,  vn_lf  ! 1st order filtered arrays          
     72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n  ! baroclinic scale factors 
     73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 
     74   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t               ! restoring period for scale factors 
     75   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv               ! restoring period for low freq. divergence 
     76   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: tildemask                 ! mask tilde tendency 
     77   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: hsm, dsm                  !  
     78   INTEGER         , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: i_int_bot 
    6779 
    6880   !! * Substitutions 
     
    8597         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   & 
    8698            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   & 
     99            &      i_int_bot(jpi,jpj), tildemask(jpi,jpj), hsm(jpi,jpj), dsm(jpi,jpj), & 
    87100            &      STAT = dom_vvl_alloc        ) 
    88101         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
     
    92105      ENDIF 
    93106      IF( ln_vvl_ztilde ) THEN 
    94          ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 
     107         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj),  & 
     108            &        hdivn_lf(jpi,jpj,jpk),  & 
     109            &           un_lf(jpi,jpj,jpk),  & 
     110            &           vn_lf(jpi,jpj,jpk),  STAT= dom_vvl_alloc ) 
    95111         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
    96112         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 
     
    124140      USE phycst,  ONLY : rpi, rsmall, rad 
    125141      !! * Local declarations 
    126       INTEGER ::   ji,jj,jk 
     142      INTEGER ::   ji, jj, jk 
    127143      INTEGER ::   ii0, ii1, ij0, ij1 
    128       REAL(wp)::   zcoef 
     144      REAL(wp)::   zcoef, zwgt, ztmp, zhmin, zhmax 
    129145      !!---------------------------------------------------------------------- 
    130146      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init') 
     
    155171      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
    156172      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
    157       CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 
     173      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' ) 
    158174      ! Vertical scale factor interpolations 
    159175      ! ------------------------------------ 
     
    199215         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    200216      END DO 
    201       hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 
    202       hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 
     217      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 
     218      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 
    203219 
    204220      ! Restoring frequencies for z_tilde coordinate 
    205221      ! ============================================ 
     222      tildemask(:,:) = 1._wp 
     223 
    206224      IF( ln_vvl_ztilde ) THEN 
    207225         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 
    208226         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp ) 
    209227         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
     228         ! tendency mask: 
     229         ! 
    210230         IF( ln_vvl_ztilde_as_zstar ) THEN 
    211231            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 
    212232            frq_rst_e3t(:,:) = 0.0_wp  
    213233            frq_rst_hdv(:,:) = 1.0_wp / rdt 
     234            tildemask(:,:) = 0._wp 
    214235         ENDIF 
     236          
    215237         IF ( ln_vvl_zstar_at_eqtor ) THEN 
    216238            DO jj = 1, jpj 
     
    219241                     ! values outside the equatorial band and transition zone (ztilde) 
    220242                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    221                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     243!                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     244               
    222245                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 
    223246                     ! values inside the equatorial band (ztilde as zstar) 
    224247                     frq_rst_e3t(ji,jj) =  0.0_wp 
    225                      frq_rst_hdv(ji,jj) =  1.0_wp / rdt 
     248!                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt 
     249                     tildemask(ji,jj) = 0._wp 
    226250                  ELSE 
    227251                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values) 
     
    229253                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    230254                        &                                          * 180._wp / 3.5_wp ) ) 
    231                      frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                & 
    232                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   & 
    233                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    234                         &                                          * 180._wp / 3.5_wp ) ) 
     255!                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                & 
     256!                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   & 
     257!                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     258!                        &                                          * 180._wp / 3.5_wp ) ) 
     259                     tildemask(ji,jj) = 0.5_wp * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     260                        &                                                 * 180._wp / 3.5_wp ) ) 
    235261                  ENDIF 
    236262               END DO 
     
    240266               ij0 = 128   ;   ij1 = 135   ;    
    241267               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp 
    242                frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt 
     268               tildemask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )   = 0.0_wp 
     269!               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt 
    243270            ENDIF 
    244271         ENDIF 
     272         ! 
     273         IF ( ln_vvl_zstar_on_shelf ) THEN 
     274            zhmin =  50._wp 
     275            zhmax = 100._wp 
     276            DO jj = 1, jpj 
     277               DO ji = 1, jpi 
     278                  zwgt = 1._wp 
     279                  IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN 
     280                     zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin) 
     281                  ELSEIF ( ht_0(ji,jj)<=zhmin) THEN 
     282                     zwgt = 0._wp 
     283                  ENDIF 
     284                  frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt) 
     285                  tildemask(ji,jj)   = MIN(tildemask(ji,jj), zwgt) 
     286               END DO 
     287            END DO 
     288         ENDIF 
     289         ! 
     290         ztmp = MAXVAL( frq_rst_hdv(:,:) ) 
     291         IF( lk_mpp )   CALL mpp_max( ztmp )                 ! max over the global domain 
     292         ! 
     293         IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' ) 
     294         ! 
    245295      ENDIF 
    246296 
     
    272322      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling. 
    273323      !!---------------------------------------------------------------------- 
    274       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 
     324      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t, ztu, ztv 
    275325      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv 
    276326      !! * Arguments 
     
    279329      !! * Local declarations 
    280330      INTEGER                                :: ji, jj, jk            ! dummy loop indices 
     331      INTEGER                                :: ib, ib_bdy, ip, jp    !   "     "     " 
    281332      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers 
    282333      REAL(wp)                               :: z2dt                  ! temporary scalars 
    283334      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars 
     335      REAL(wp)                               :: zalpha, zwgt          ! temporary scalars 
     336      REAL(wp)                               :: zdu, zdv 
    284337      LOGICAL                                :: ll_do_bclinic         ! temporary logical 
    285338      !!---------------------------------------------------------------------- 
    286339      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt') 
    287340      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 
    288       CALL wrk_alloc( jpi, jpj, jpk, ze3t                    ) 
     341      CALL wrk_alloc( jpi, jpj, jpk, ze3t, ztu, ztv ) 
    289342 
    290343      IF(kt == nit000)   THEN 
     
    296349      ll_do_bclinic = .TRUE. 
    297350      IF( PRESENT(kcall) ) THEN 
    298          IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 
     351         IF ( kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE. 
    299352      ENDIF 
    300353 
     
    307360      !                                                ! --------------------------------------------- ! 
    308361 
    309       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     362      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - ssmask(:,:) ) 
    310363      DO jk = 1, jpkm1 
    311364         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) 
    312365         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    313366      END DO 
    314  
    315       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    316          !                                                            ! ------baroclinic part------ ! 
    317  
    318          ! I - initialization 
    319          ! ================== 
    320  
    321          ! 1 - barotropic divergence 
    322          ! ------------------------- 
     367       
     368      IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     369          
     370         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms 
     371         un_td(:,:,:) = 0.0_wp        ! Transport corrections 
     372         vn_td(:,:,:) = 0.0_wp 
     373 
    323374         zhdiv(:,:) = 0. 
    324          zht(:,:)   = 0. 
    325375         DO jk = 1, jpkm1 
    326376            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 
    327             zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
    328          END DO 
    329          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
    330  
    331          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only) 
    332          ! -------------------------------------------------- 
    333          IF( ln_vvl_ztilde ) THEN 
    334             IF( kt .GT. nit000 ) THEN 
     377         END DO 
     378         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )       
     379 
     380         ! Thickness advection: 
     381         ! -------------------- 
     382         ! Set advection velocities and source term 
     383         IF ( ln_vvl_ztilde ) THEN 
     384            ! 
     385            IF ((kt==nit000).AND.(neuler==0)) THEN 
    335386               DO jk = 1, jpkm1 
    336                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
    337                      &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
     387                  ztu(:,:,jk) = un(:,:,jk) 
     388                  ztv(:,:,jk) = vn(:,:,jk) 
     389               END DO 
     390            ELSE 
     391               DO jk = 1, jpkm1 
     392                  ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk) 
     393                  ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk) 
     394                  tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk) 
    338395               END DO 
    339396            ENDIF 
    340          END IF 
    341  
    342          ! II - after z_tilde increments of vertical scale factors 
    343          ! ======================================================= 
    344          tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms 
    345  
    346          ! 1 - High frequency divergence term 
    347          ! ---------------------------------- 
    348          IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
     397            ! 
     398         ELSEIF ( ln_vvl_layer ) THEN 
     399            ! 
    349400            DO jk = 1, jpkm1 
    350                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    351             END DO 
    352          ELSE                         ! layer case 
     401               ztu(:,:,jk) = un(:,:,jk) 
     402               ztv(:,:,jk) = vn(:,:,jk) 
     403            END DO 
     404            ! 
     405         ENDIF 
     406         ! 
     407         ! Do advection 
     408         IF     (ln_vvl_adv_fct) THEN 
     409            CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv ) 
     410            ! 
     411         ELSEIF (ln_vvl_adv_cn2) THEN 
    353412            DO jk = 1, jpkm1 
    354                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    355             END DO 
    356          END IF 
    357  
    358          ! 2 - Restoring term (z-tilde case only) 
    359          ! ------------------ 
     413               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * hdivn(:,:,jk) 
     414            END DO 
     415         ENDIF 
     416         !  
     417         ! Thickness anomlaly diffusion: 
     418         ! ----------------------------- 
     419         zwu(:,:) = 0.0_wp 
     420         zwv(:,:) = 0.0_wp 
     421         ztu(:,:,:) = 0.0_wp 
     422         ztv(:,:,:) = 0.0_wp 
     423 
     424         IF ( ln_vvl_blp ) THEN  ! Bilaplacian 
     425            DO jk = 1, jpkm1 
     426               DO jj = 1, jpjm1                 ! First derivative (gradient) 
     427                  DO ji = 1, fs_jpim1   ! vector opt.                   
     428                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     429                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     430                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     431                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     432                  END DO 
     433               END DO 
     434 
     435               DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
     436                  DO ji = fs_2, fs_jpim1 ! vector opt. 
     437                     zht(ji,jj) =  rn_ahe3_blp * r1_e12t(ji,jj) * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
     438                        &                                           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   ) 
     439                  END DO 
     440               END DO 
     441 
     442#if defined key_bdy 
     443               DO ib_bdy=1, nb_bdy 
     444                  DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 
     445                     ji = idx_bdy(ib_bdy)%nbi(ib,1) 
     446                     jj = idx_bdy(ib_bdy)%nbj(ib,1) 
     447                     zht(ji,jj) = 0._wp 
     448                  END DO  
     449               END DO 
     450#endif 
     451               CALL lbc_lnk( zht, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn) 
     452 
     453               DO jj = 1, jpjm1                 ! third derivative (gradient) 
     454                  DO ji = 1, fs_jpim1   ! vector opt. 
     455                     ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) ) 
     456                     ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( zht(ji  ,jj+1) - zht(ji,jj) ) 
     457                     zwu(ji,jj) = zwu(ji,jj) + ztu(ji,jj,jk) 
     458                     zwv(ji,jj) = zwv(ji,jj) + ztv(ji,jj,jk) 
     459                  END DO 
     460               END DO 
     461            END DO 
     462         ENDIF 
     463 
     464         IF ( ln_vvl_lap ) THEN    ! Laplacian 
     465            DO jk = 1, jpkm1                    ! First derivative (gradient) 
     466               DO jj = 1, jpjm1 
     467                  DO ji = 1, fs_jpim1   ! vector opt.                   
     468                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     469                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     470                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     471                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     472                     zwu(ji,jj) = zwu(ji,jj) + zdu 
     473                     zwv(ji,jj) = zwv(ji,jj) + zdv 
     474                     ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu 
     475                     ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv 
     476                  END DO 
     477               END DO 
     478            END DO 
     479         ENDIF  
     480 
     481         ! Ensure barotropic fluxes are null: 
     482!         DO jj = 1, jpj 
     483!            DO ji = 1, jpi 
     484!               DO jk = 1, jpkm1 
     485!                  ztu(ji,jj,jk) = ztu(ji,jj,jk) - zwu(ji,jj)*fse3u_b(ji,jj,jk)*hur_b(ji,jj)*umask(ji,jj,jk) 
     486!                  ztv(ji,jj,jk) = ztv(ji,jj,jk) - zwv(ji,jj)*fse3v_b(ji,jj,jk)*hvr_b(ji,jj)*vmask(ji,jj,jk) 
     487!               END DO 
     488!            END DO 
     489!         END DO 
     490         DO jj = 1, jpj 
     491            DO ji = 1, jpi 
     492               ztu(ji,jj,mbku(ji,jj)) = ztu(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
     493               ztv(ji,jj,mbkv(ji,jj)) = ztv(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
     494            END DO 
     495         END DO 
     496 
     497         ! divergence of diffusive fluxes 
     498         DO jk = 1, jpkm1 
     499            DO jj = 2, jpjm1 
     500               DO ji = fs_2, fs_jpim1   ! vector opt. 
     501                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk) - ztu(ji,jj,jk)    & 
     502                     &                                          +     ztv(ji  ,jj-1,jk) - ztv(ji,jj,jk)    & 
     503                     &                                            ) * r1_e12t(ji,jj) 
     504               END DO 
     505            END DO 
     506         END DO 
     507  
     508         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:) 
     509         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:) 
     510         CALL lbc_lnk( un_td , 'U' , -1.) 
     511         CALL lbc_lnk( vn_td , 'V' , -1.)   
     512         ! 
     513         ! 
     514         ! Restoring: 
    360515         IF( ln_vvl_ztilde ) THEN 
    361516            DO jk = 1, jpk 
    362517               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
    363518            END DO 
    364          END IF 
    365  
    366          ! 3 - Thickness diffusion term 
    367          ! ---------------------------- 
    368          zwu(:,:) = 0.0_wp 
    369          zwv(:,:) = 0.0_wp 
    370          ! a - first derivative: diffusive fluxes 
     519         ENDIF 
     520 
     521         ! Remove "external thickness" tendency: 
    371522         DO jk = 1, jpkm1 
    372             DO jj = 1, jpjm1 
    373                DO ji = 1, fs_jpim1   ! vector opt. 
    374                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
    375                                   & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    376                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
    377                                   & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    378                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    379                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    380                END DO 
    381             END DO 
    382          END DO 
    383          ! b - correction for last oceanic u-v points 
    384          DO jj = 1, jpj 
    385             DO ji = 1, jpi 
    386                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    387                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    388             END DO 
    389          END DO 
    390          ! c - second derivative: divergence of diffusive fluxes 
    391          DO jk = 1, jpkm1 
    392             DO jj = 2, jpjm1 
    393                DO ji = fs_2, fs_jpim1   ! vector opt. 
    394                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    395                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    396                      &                                            ) * r1_e12t(ji,jj) 
    397                END DO 
    398             END DO 
    399          END DO 
    400          ! d - thickness diffusion transport: boundary conditions 
    401          !     (stored for tracer advction and continuity equation) 
    402          CALL lbc_lnk( un_td , 'U' , -1._wp) 
    403          CALL lbc_lnk( vn_td , 'V' , -1._wp) 
    404  
    405          ! 4 - Time stepping of baroclinic scale factors 
    406          ! --------------------------------------------- 
     523            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) +   fse3t_n(:,:,jk) * zhdiv(:,:) 
     524         END DO  
     525                    
    407526         ! Leapfrog time stepping 
    408527         ! ~~~~~~~~~~~~~~~~~~~~~~ 
     
    412531            z2dt = 2.0_wp * rdt 
    413532         ENDIF 
    414          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 
     533 
    415534         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
     535 
     536         ! Revert to zstar locally: 
     537         ! ~~~~~~~~~~~~~~~~~~~~~~~~ 
     538         DO jk=1,jpkm1 
     539            tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) * tildemask(:,:)  
     540         END DO 
     541 
     542         ! Ensure layer separation: 
     543         ! ~~~~~~~~~~~~~~~~~~~~~~~~ 
     544         IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt ) 
     545   
     546         ! Boundary conditions: 
     547         ! ~~~~~~~~~~~~~~~~~~~~ 
     548#if defined key_bdy 
     549         DO ib_bdy=1, nb_bdy 
     550            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 
     551!!            DO ib = 1, idx_bdy(ib_bdy)%nblen(1) 
     552               ji = idx_bdy(ib_bdy)%nbi(ib,1) 
     553               jj = idx_bdy(ib_bdy)%nbj(ib,1) 
     554               zwgt = idx_bdy(ib_bdy)%nbw(ib,1) 
     555               ip = bdytmask(ji+1,jj  ) - bdytmask(ji-1,jj  ) 
     556               jp = bdytmask(ji  ,jj+1) - bdytmask(ji  ,jj-1) 
     557               DO jk = 1, jpkm1   
     558                  tilde_e3t_a(ji,jj,jk) = 0.e0 
     559!!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt) 
     560!!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk) 
     561               END DO 
     562            END DO  
     563         END DO 
     564#endif 
     565         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )         
    416566 
    417567         ! Maximum deformation control 
     
    446596            ENDIF 
    447597         ENDIF 
    448          ! - ML - end test 
    449          ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 
    450          tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) ) 
    451          tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 
    452  
    453          ! 
    454          ! "tilda" change in the after scale factor 
    455          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
     598 
     599         IF ( ln_vvl_ztilde ) THEN 
     600            zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
     601            DO jk = 1, jpkm1 
     602               ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:) + un_td(:,:,jk) 
     603               ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk) 
     604               ze3t(:,:,jk) = -fse3t_n(:,:,jk) * zhdiv(:,:) 
     605! 
     606!               DO jj = 2, jpjm1 
     607!                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     608!                      
     609!                     ze3t(ji,jj,jk) =   -fse3t_n(ji,jj,jk) * zhdiv(ji,jj)   & 
     610!                        & + (  un_td(ji,jj,jk) - un_td(ji-1,jj  ,jk)        & 
     611!                        &    + vn_td(ji,jj,jk) - vn_td(ji  ,jj-1,jk)  )     & 
     612!                        & / ( e1t(ji,jj) * e2t(ji,jj) ) 
     613!                  END DO 
     614!               END DO  
     615            END DO 
     616            ! 
     617            un_lf(:,:,:)    = un_lf(:,:,:)    * (1._wp - zalpha) + zalpha * ztu(:,:,:) 
     618            vn_lf(:,:,:)    = vn_lf(:,:,:)    * (1._wp - zalpha) + zalpha * ztv(:,:,:) 
     619            hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:) 
     620         ENDIF 
     621 
     622      ENDIF 
     623 
     624      IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.(.NOT.ll_do_bclinic) ) THEN 
     625         zhdiv(:,:) = 0. 
    456626         DO jk = 1, jpkm1 
    457             dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 
    458          END DO 
    459          ! III - Barotropic repartition of the sea surface height over the baroclinic profile 
    460          ! ================================================================================== 
    461          ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 
    462          ! - ML - baroclinicity error should be better treated in the future 
    463          !        i.e. locally and not spread over the water column. 
    464          !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 
    465          zht(:,:) = 0. 
     627            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * (hdivn(:,:,jk) - hdivb(:,:,jk)) 
     628         END DO 
     629         zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
     630 
     631         IF( neuler == 0 .AND. kt == nit000 ) THEN 
     632            z2dt =  rdt 
     633         ELSE 
     634            z2dt = 2.0_wp * rdt 
     635         ENDIF 
     636 
    466637         DO jk = 1, jpkm1 
    467             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    468          END DO 
    469          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     638            tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * &  
     639                                  & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:)) 
     640         END DO 
     641      ENDIF 
     642 
     643      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
     644      !                                             ! ---baroclinic part--------- ! 
    470645         DO jk = 1, jpkm1 
    471             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    472          END DO 
    473  
     646            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                      
     647         END DO 
    474648      ENDIF 
    475649 
    476       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
    477       !                                           ! ---baroclinic part--------- ! 
     650      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging 
     651         ! 
     652         zht(:,:) = 0.0_wp 
    478653         DO jk = 1, jpkm1 
    479             fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    480          END DO 
    481       ENDIF 
    482  
    483       IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging 
    484          ! 
     654            zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
     655         END DO 
    485656         IF( lwp ) WRITE(numout, *) 'kt =', kt 
    486657         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    487             z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 
     658            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 ) 
    488659            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain 
    489660            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 
    490661         END IF 
    491662         ! 
     663         z_tmin = MINVAL( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0  ) 
     664         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
     665         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin          
     666         ! 
     667         z_tmin = MINVAL( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 ) 
     668         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
     669         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3u_n) =', z_tmin 
     670         ! 
     671         z_tmin = MINVAL( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 ) 
     672         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
     673         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3v_n) =', z_tmin 
     674         ! 
     675         z_tmin = MINVAL( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 ) 
     676         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
     677         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3f_n) =', z_tmin 
     678         ! 
    492679         zht(:,:) = 0.0_wp 
    493680         DO jk = 1, jpkm1 
     
    496683         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 
    497684         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    498          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax 
     685         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn  -SUM(fse3t_n))) =', z_tmax 
    499686         ! 
    500687         zht(:,:) = 0.0_wp 
    501688         DO jk = 1, jpkm1 
    502             zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk) 
    503          END DO 
    504          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 
     689            zht(:,:) = zht(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
     690         END DO 
     691         zwu(:,:) = 0._wp 
     692         DO jj = 1, jpjm1 
     693            DO ji = 1, fs_jpim1   ! vector opt. 
     694               zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e12u(ji,jj)                             & 
     695                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji+1,jj) * sshn(ji+1,jj) ) 
     696            END DO 
     697         END DO 
     698         CALL lbc_lnk( zwu(:,:), 'U', 1._wp ) 
     699         z_tmax = MAXVAL( umask(:,:,1) * umask_i(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) ) 
    505700         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    506          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax 
     701         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(fse3u_n))) =', z_tmax 
    507702         ! 
    508703         zht(:,:) = 0.0_wp 
    509704         DO jk = 1, jpkm1 
    510             zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk) 
    511          END DO 
    512          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 
     705            zht(:,:) = zht(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
     706         END DO 
     707         zwv(:,:) = 0._wp 
     708         DO jj = 1, jpjm1 
     709            DO ji = 1, fs_jpim1   ! vector opt. 
     710               zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e12v(ji,jj)                             & 
     711                        &          * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji,jj+1) * sshn(ji,jj+1) ) 
     712            END DO 
     713         END DO 
     714         CALL lbc_lnk( zwv(:,:), 'V', 1._wp ) 
     715         z_tmax = MAXVAL( vmask(:,:,1) * vmask_i(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) ) 
    513716         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    514          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax 
    515          ! 
    516          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) ) 
     717         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(fse3v_n))) =', z_tmax 
     718         ! 
     719         zht(:,:) = 0.0_wp 
     720         DO jk = 1, jpkm1 
     721            DO jj = 1, jpjm1 
     722               zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk) 
     723            END DO 
     724         END DO 
     725         zwu(:,:) = 0._wp 
     726         DO jj = 1, jpjm1 
     727            DO ji = 1, fs_jpim1   ! vector opt. 
     728               zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e12f(ji,jj)  & 
     729                        &          * (  e12t(ji  ,jj) * sshn(ji  ,jj) + e12t(ji  ,jj+1) * sshn(ji  ,jj+1) & 
     730                        &             + e12t(ji+1,jj) * sshn(ji+1,jj) + e12t(ji+1,jj+1) * sshn(ji+1,jj+1) ) 
     731            END DO 
     732         END DO 
     733         CALL lbc_lnk( zht(:,:), 'F', 1._wp ) 
     734         CALL lbc_lnk( zwu(:,:), 'F', 1._wp ) 
     735         z_tmax = MAXVAL( fmask(:,:,1) * fmask_i(:,:) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) ) 
    517736         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    518          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 
    519          ! 
    520          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) ) 
    521          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    522          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 
    523          ! 
    524          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) ) 
    525          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    526          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 
     737         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(fse3f_n))) =', z_tmax 
     738         ! 
    527739      END IF 
    528740 
     
    549761 
    550762      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 
    551       CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     ) 
     763      CALL wrk_dealloc( jpi, jpj, jpk, ze3t, ztu, ztv           ) 
    552764 
    553765      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt') 
     
    583795      INTEGER, INTENT( in )               :: kt       ! time step 
    584796      !! * Local declarations 
    585       INTEGER                             :: ji,jj,jk       ! dummy loop indices 
    586       REAL(wp)                            :: zcoef 
     797      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwrk 
     798      INTEGER                             :: jk       ! dummy loop indices 
    587799      !!---------------------------------------------------------------------- 
    588800 
     
    607819         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    608820      ENDIF 
     821 
    609822      fsdept_b(:,:,:) = fsdept_n(:,:,:) 
    610823      fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 
     
    620833      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt 
    621834      ! - JC - hu_b, hv_b, hur_b, hvr_b also 
    622       CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  ) 
     835      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n (:,:,:), 'F'  ) 
    623836      ! Vertical scale factor interpolations 
    624837      ! ------------------------------------ 
     
    631844      ! t- and w- points depth 
    632845      ! ---------------------- 
    633       ! set the isf depth as it is in the initial step 
    634846      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    635847      fsdepw_n(:,:,1) = 0.0_wp 
    636848      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
    637  
    638849      DO jk = 2, jpk 
    639          DO jj = 1,jpj 
    640             DO ji = 1,jpi 
    641               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    642                                                                  ! 1 for jk = mikt 
    643                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    644                fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    645                fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
    646                    &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk))  
    647                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 
    648             END DO 
    649          END DO 
     850         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 
     851         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 
     852         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:) 
    650853      END DO 
    651  
    652854      ! Local depth and Inverse of the local depth of the water column at u- and v- points 
    653855      ! ---------------------------------------------------------------------------------- 
     
    666868      END DO 
    667869 
     870      ! Write additional diagnostics  
     871      ! ============================ 
     872      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) CALL dom_vvl_dia( kt)  
     873 
    668874      ! write restart file 
    669875      ! ================== 
     
    674880   END SUBROUTINE dom_vvl_sf_swp 
    675881 
     882   SUBROUTINE dom_vvl_dia( kt ) 
     883      !!---------------------------------------------------------------------- 
     884      !!                ***  ROUTINE dom_vvl_dia  *** 
     885      !!                    
     886      !! ** Purpose :  Output some diagnostics in ztilde/zlayer cases 
     887      !! 
     888      !!---------------------------------------------------------------------- 
     889      !! * Arguments 
     890      INTEGER, INTENT( in )               :: kt       ! time step 
     891      !! * Local declarations 
     892      INTEGER                             :: ji,jj,jk ! dummy loop indices 
     893      REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, ztmp1, z2dt 
     894      REAL(wp), DIMENSION(4)              :: zr1 
     895      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zout 
     896      !!---------------------------------------------------------------------- 
     897      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_dia') 
     898      ! 
     899      CALL wrk_alloc( jpi, jpj, jpk, zwdw, zout ) 
     900      ! 
     901      ! Compute internal interfaces depths: 
     902      !------------------------------------ 
     903      IF ( iom_use("dh_tilde").OR.iom_use("depw_tilde").OR.iom_use("stiff_tilde")) THEN 
     904         zwdw(:,:,1) = 0.e0 
     905         DO jj = 1, jpj 
     906            DO ji = 1, jpi 
     907               DO jk = 2, jpkm1 
     908                  zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + &  
     909                                 & (tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     910               END DO 
     911            END DO 
     912         END DO 
     913      ENDIF 
     914      ! 
     915      ! Output interface depth anomaly: 
     916      ! ------------------------------- 
     917      IF ( iom_use("depw_tilde") ) CALL iom_put( "depw_tilde", (zwdw(:,:,:)-gdepw_0(:,:,:))*tmask(:,:,:) ) 
     918      ! 
     919      ! Output grid stiffness (T-points): 
     920      ! --------------------------------- 
     921      IF ( iom_use("stiff_tilde"  ) ) THEN 
     922         zr1(:)   = 0.e0 
     923         zout(:,:,jpk) = 0.e0       
     924         DO ji = 2, jpim1 
     925            DO jj = 2, jpjm1 
     926               DO jk = 1, jpkm1 
     927                  zr1(1) = umask(ji-1,jj  ,jk) *abs( (zwdw(ji  ,jj  ,jk  )-zwdw(ji-1,jj  ,jk  )  &  
     928                       &                             +zwdw(ji  ,jj  ,jk+1)-zwdw(ji-1,jj  ,jk+1)) & 
     929                       &                            /(zwdw(ji  ,jj  ,jk  )+zwdw(ji-1,jj  ,jk  )  & 
     930                       &                             -zwdw(ji  ,jj  ,jk+1)-zwdw(ji-1,jj  ,jk+1) + rsmall) ) 
     931                  zr1(2) = umask(ji  ,jj  ,jk) *abs( (zwdw(ji+1,jj  ,jk  )-zwdw(ji  ,jj  ,jk  )  & 
     932                       &                             +zwdw(ji+1,jj  ,jk+1)-zwdw(ji  ,jj  ,jk+1)) & 
     933                       &                            /(zwdw(ji+1,jj  ,jk  )+zwdw(ji  ,jj  ,jk  )  & 
     934                       &                             -zwdw(ji+1,jj  ,jk+1)-zwdw(ji  ,jj  ,jk+1) + rsmall) ) 
     935                  zr1(3) = vmask(ji  ,jj  ,jk) *abs( (zwdw(ji  ,jj+1,jk  )-zwdw(ji  ,jj  ,jk  )  & 
     936                       &                             +zwdw(ji  ,jj+1,jk+1)-zwdw(ji  ,jj  ,jk+1)) & 
     937                       &                            /(zwdw(ji  ,jj+1,jk  )+zwdw(ji  ,jj  ,jk  )  & 
     938                       &                             -zwdw(ji  ,jj+1,jk+1)-zwdw(ji  ,jj  ,jk+1) + rsmall) ) 
     939                  zr1(4) = vmask(ji  ,jj-1,jk) *abs( (zwdw(ji  ,jj  ,jk  )-zwdw(ji  ,jj-1,jk  )  & 
     940                       &                             +zwdw(ji  ,jj  ,jk+1)-zwdw(ji  ,jj-1,jk+1)) & 
     941                       &                            /(zwdw(ji  ,jj  ,jk  )+zwdw(ji  ,jj-1,jk  )  & 
     942                       &                             -zwdw(ji,  jj  ,jk+1)-zwdw(ji  ,jj-1,jk+1) + rsmall) ) 
     943                  zout(ji,jj,jk) = MAXVAL(zr1(1:4)) 
     944               END DO 
     945            END DO 
     946         END DO 
     947 
     948         CALL lbc_lnk( zout, 'T', 1. ) 
     949         CALL iom_put( "stiff_tilde", zout(:,:,:) )  
     950      END IF 
     951      ! Output Horizontal Laplacian of interfaces depths (W-points): 
     952      ! ------------------------------------------------------------ 
     953      IF ( iom_use("dh_tilde")   ) THEN 
     954         ! 
     955         zout(:,:,1  )=0._wp 
     956         zout(:,:,jpk)=0._wp 
     957         DO jk = 2, jpkm1 
     958            DO jj = 2, jpjm1 
     959               DO ji = fs_2, fs_jpim1   ! vector opt. 
     960                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) ) 
     961                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) ) 
     962                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) ) 
     963                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) ) 
     964!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj))   
     965                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj)  
     966!                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk) 
     967                  zout(ji,jj,jk) = ztmp1 
     968               END DO 
     969            END DO             
     970         END DO 
     971         ! Mask open boundaries: 
     972#if defined key_bdy 
     973         IF (lk_bdy) THEN 
     974            DO jk = 1, jpkm1 
     975               zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:) 
     976            END DO 
     977         ENDIF 
     978#endif 
     979         CALL lbc_lnk( zout(:,:,:), 'T', 1. ) 
     980         zwdw(:,:,:) = zout(:,:,:) 
     981         DO jk = 2, jpkm1 
     982            DO jj = 2, jpjm1 
     983               DO ji = fs_2, fs_jpim1   ! vector opt. 
     984                  zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji  ,jj  ,jk) ) 
     985                  zufi   = umask(ji  ,jj,jk) * re2u_e1u(ji  ,jj) * ( zwdw(ji  ,jj,jk) - zwdw(ji+1,jj  ,jk) ) 
     986                  zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji  ,jj  ,jk) ) 
     987                  zvfj   = vmask(ji,jj  ,jk) * re1v_e2v(ji,jj  ) * ( zwdw(ji,jj  ,jk) - zwdw(ji  ,jj+1,jk) ) 
     988!                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj))   
     989                  ztmp1  = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj)  
     990                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk) 
     991               END DO 
     992            END DO             
     993         END DO 
     994         ! Mask open boundaries: 
     995#if defined key_bdy 
     996         IF (lk_bdy) THEN 
     997            DO jk = 1, jpkm1 
     998               zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:) 
     999            END DO 
     1000         ENDIF 
     1001#endif 
     1002         CALL lbc_lnk( zout(:,:,:), 'T', 1. ) 
     1003         ! 
     1004         CALL iom_put( "dh_tilde", zout(:,:,:) ) 
     1005      ENDIF 
     1006      ! 
     1007      ! Output vertical Laplacian of interfaces depths (W-points): 
     1008      ! ---------------------------------------------------------- 
     1009      IF ( iom_use("dz_tilde"  ) ) THEN 
     1010         zout(:,:,1  ) = 0._wp 
     1011         zout(:,:,jpk) = 0._wp 
     1012         DO jk=2,jpkm1 
     1013            zout(:,:,jk) = 2._wp*ABS(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)-tilde_e3t_n(:,:,jk-1)-e3t_0(:,:,jk-1)) &  
     1014                               &   /(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)+tilde_e3t_n(:,:,jk-1)+e3t_0(:,:,jk-1)) & 
     1015                               &   * tmask(:,:,jk) 
     1016         END DO 
     1017         CALL iom_put( "dz_tilde", zout(:,:,:) )  
     1018      END IF 
     1019      ! 
     1020      ! 
     1021      ! Output low pass U-velocity: 
     1022      ! --------------------------- 
     1023      IF ( iom_use("un_lf_tilde"  ).AND.ln_vvl_ztilde ) THEN 
     1024         zout(:,:,jpk) = 0.e0   
     1025         DO jk=1,jpkm1 
     1026            zout(:,:,jk) = un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:) 
     1027         END DO 
     1028         CALL iom_put( "un_lf_tilde", zout(:,:,:) ) 
     1029      END IF 
     1030      ! 
     1031      ! Output low pass V-velocity: 
     1032      ! --------------------------- 
     1033      IF ( iom_use("vn_lf_tilde"  ).AND.ln_vvl_ztilde ) THEN 
     1034         zout(:,:,jpk) = 0.e0   
     1035         DO jk=1,jpkm1 
     1036            zout(:,:,jk) = vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:) 
     1037         END DO 
     1038         CALL iom_put( "vn_lf_tilde", zout(:,:,:) ) 
     1039      END IF 
     1040      ! 
     1041      CALL wrk_dealloc( jpi, jpj, jpk, zwdw, zout ) 
     1042      ! 
     1043      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_dia') 
     1044      ! 
     1045   END SUBROUTINE dom_vvl_dia 
    6761046 
    6771047   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 
     
    6911061      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    6921062      !! * Local declarations 
     1063      INTEGER ::   nmet                                                ! horizontal interpolation method 
    6931064      INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
     1065      INTEGER ::   jkbot                                               ! bottom level 
    6941066      LOGICAL ::   l_is_orca                                           ! local logical 
     1067      REAL(wp) :: zmin, zdo, zup, ztap 
     1068      REAL(wp), POINTER, DIMENSION(:,:)   :: zs                        ! Surface interface depth anomaly 
     1069      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw                        ! Interface depth anomaly 
    6951070      !!---------------------------------------------------------------------- 
    6961071      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol') 
     
    6981073      l_is_orca = .FALSE. 
    6991074      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations 
     1075 
     1076  
     1077      nmet=1 ! Original method (Surely wrong) 
     1078!      nmet=1 ! Interface interpolation 
     1079!      nmet=2 ! Internal interfaces interpolation only, spread barotropic increment 
     1080!      nmet=1 
     1081      ztap = 0.1_wp ! Minimum fraction of T-point thickness at cell interfaces 
     1082 
     1083      IF ( (nmet==1).OR.(nmet==2) ) THEN 
     1084         SELECT CASE ( pout ) 
     1085            ! 
     1086         CASE( 'U', 'V', 'F' ) 
     1087            ! Compute interface depth anomaly at T-points 
     1088            CALL wrk_alloc( jpi, jpj, jpk, zw ) 
     1089            CALL wrk_alloc( jpi, jpj, zs ) 
     1090            ! 
     1091            zw(:,:,:) =  0._wp     
     1092            ! 
     1093            DO jj = 1, jpj 
     1094               DO ji = 1, jpi 
     1095                  jkbot = mbkt(ji,jj) 
     1096                  DO jk=jkbot,1,-1 
     1097                     zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk) 
     1098                  END DO 
     1099               END DO 
     1100            END DO  
     1101            ! 
     1102            ! 
     1103            IF (nmet==2) THEN        ! Consider "internal" interfaces only 
     1104               zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh) 
     1105               ! 
     1106               zw(:,:,:) = 0._wp 
     1107               DO jj = 1, jpj 
     1108                  DO ji = 1, jpi 
     1109                     jkbot = mbkt(ji,jj) 
     1110                     DO jk=jkbot,1,-1 
     1111                        zw(ji,jj,jk) = zw(ji,jj,jk+1) - ( pe3_in(ji,jj,jk)                     & 
     1112                                     & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1))  & 
     1113                                     & - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
     1114                     END DO 
     1115                  END DO 
     1116               END DO  
     1117            ENDIF  
     1118            ! 
     1119         END SELECT 
     1120      END IF 
     1121 
     1122      pe3_out(:,:,:) = 0._wp 
    7001123 
    7011124      SELECT CASE ( pout ) 
    7021125         !               ! ------------------------------------- ! 
    7031126      CASE( 'U' )        ! interpolation from T-point to U-point ! 
    704          !               ! ------------------------------------- ! 
     1127         !               ! ------------------------------------- !      
    7051128         ! horizontal surface weighted interpolation 
    706          DO jk = 1, jpk 
     1129         IF (nmet==0) THEN 
     1130            DO jk = 1, jpk 
     1131               DO jj = 1, jpjm1 
     1132                  DO ji = 1, fs_jpim1   ! vector opt. 
     1133                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
     1134                        &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     1135                        &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     1136                  END DO 
     1137               END DO 
     1138            END DO 
     1139         ENDIF 
     1140         ! 
     1141         IF ( (nmet==1).OR.(nmet==2) ) THEN 
    7071142            DO jj = 1, jpjm1 
    7081143               DO ji = 1, fs_jpim1   ! vector opt. 
    709                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
    710                      &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    711                      &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    712                END DO 
    713             END DO 
    714          END DO 
     1144                  ! Correction at last level: 
     1145                  jkbot = mbku(ji,jj) 
     1146                  pe3_out(ji,jj,jkbot) = -0.5_wp * umask(ji,jj,jkbot) * r1_e12u(ji,jj) & 
     1147                        &                    * (   e12t(ji  ,jj) * zw(ji  ,jj,jkbot)   & 
     1148                        &                        + e12t(ji+1,jj) * zw(ji+1,jj,jkbot)   ) 
     1149                  ! 
     1150                  ! If there is a step, taper bottom interface: 
     1151                  IF (hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ) THEN 
     1152                     IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN 
     1153!                        zmin = ztap * pe3_in(ji+1,jj,jkbot) 
     1154                        zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot)) 
     1155                     ELSE 
     1156!                        zmin = ztap * pe3_in(ji  ,jj,jkbot) 
     1157                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot)) 
     1158                     ENDIF 
     1159                     zmin = -e3u_0(ji,jj,jkbot) + zmin 
     1160                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin) 
     1161                  ENDIF 
     1162                  ! 
     1163                  zdo =  -pe3_out(ji,jj,jkbot) 
     1164                  DO jk=jkbot-1,1,-1 
     1165                     zup = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji  ,jj) & 
     1166                              &                     *(   e12t(ji  ,jj) * zw(ji  ,jj,jk)  & 
     1167                              &                         +e12t(ji+1,jj) * zw(ji+1,jj,jk)  ) 
     1168                     pe3_out(ji,jj,jk) = zdo - zup   
     1169                     zdo =  zdo - pe3_out(ji,jj,jk) 
     1170                  END DO 
     1171               END DO 
     1172            END DO 
     1173         END IF 
     1174         !  
     1175         IF (nmet==2) THEN           ! Spread sea level anomaly 
     1176            DO jj = 1, jpjm1 
     1177               DO ji = 1, fs_jpim1   ! vector opt. 
     1178                  DO jk=1,jpk 
     1179                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                 & 
     1180                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )             &  
     1181                           &               / ( hu_0(ji,jj) + 1._wp - umask_i(ji,jj) )            & 
     1182                           &               * 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)           & 
     1183                           &               * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji+1,jj)*zs(ji+1,jj) )        
     1184                  END DO 
     1185               END DO 
     1186            END DO 
     1187            ! 
     1188         ENDIF 
    7151189         ! 
    7161190         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
     
    7181192         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 
    7191193         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
     1194         !  
     1195         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs ) 
     1196         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw ) 
     1197         ! 
    7201198         !               ! ------------------------------------- ! 
    7211199      CASE( 'V' )        ! interpolation from T-point to V-point ! 
    7221200         !               ! ------------------------------------- ! 
    7231201         ! horizontal surface weighted interpolation 
    724          DO jk = 1, jpk 
     1202         IF (nmet==0) THEN 
     1203            DO jk = 1, jpk 
     1204               DO jj = 1, jpjm1 
     1205                  DO ji = 1, fs_jpim1   ! vector opt. 
     1206                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
     1207                        &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     1208                        &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     1209                  END DO 
     1210               END DO 
     1211            END DO 
     1212         ENDIF 
     1213         ! 
     1214         IF ( (nmet==1).OR.(nmet==2) ) THEN 
    7251215            DO jj = 1, jpjm1 
    7261216               DO ji = 1, fs_jpim1   ! vector opt. 
    727                   pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
    728                      &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    729                      &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    730                END DO 
    731             END DO 
    732          END DO 
     1217                  ! Correction at last level: 
     1218                  jkbot = mbkv(ji,jj) 
     1219                  pe3_out(ji,jj,jkbot) = -0.5_wp * vmask(ji,jj,jkbot) * r1_e12v(ji,jj) & 
     1220                        &                    * (   e12t(ji,jj  ) * zw(ji,jj  ,jkbot)   & 
     1221                        &                        + e12t(ji,jj+1) * zw(ji,jj+1,jkbot)   ) 
     1222                  ! 
     1223                  ! If there is a step, taper bottom interface: 
     1224                  IF (hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ) THEN 
     1225                     IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN 
     1226!                        zmin = ztap * pe3_in(ji,jj+1,jkbot) 
     1227                        zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot)) 
     1228                     ELSE 
     1229!                        zmin = ztap * pe3_in(ji,jj  ,jkbot) 
     1230                        zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot)) 
     1231                     ENDIF 
     1232                     zmin = -e3v_0(ji,jj,jkbot) + zmin 
     1233                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin) 
     1234                  ENDIF 
     1235                  ! 
     1236                  zdo =  -pe3_out(ji,jj,jkbot) 
     1237                  DO jk=jkbot-1,1,-1 
     1238                     zup = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj  ) & 
     1239                              &                     *  ( e12t(ji,jj  ) * zw(ji,jj  ,jk) & 
     1240                              &                         +e12t(ji,jj+1) * zw(ji,jj+1,jk) ) 
     1241                     ! 
     1242                     pe3_out(ji,jj,jk) = zdo - zup  
     1243                     zdo =  zdo - pe3_out(ji,jj,jk) 
     1244                  END DO 
     1245               END DO 
     1246            END DO 
     1247         END IF 
     1248         !  
     1249         IF (nmet==2) THEN           ! Spread sea level anomaly 
     1250            ! 
     1251            DO jj = 1, jpjm1 
     1252               DO ji = 1, fs_jpim1   ! vector opt. 
     1253                  DO jk=1,jpk 
     1254                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                 & 
     1255                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )             &  
     1256                           &               / ( hv_0(ji,jj) + 1._wp - vmask_i(ji,jj) )            & 
     1257                           &               * 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)           & 
     1258                           &               * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji,jj+1)*zs(ji,jj+1) )        
     1259                  END DO 
     1260               END DO 
     1261            END DO 
     1262            ! 
     1263         ENDIF 
    7331264         ! 
    7341265         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
     
    7361267         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 
    7371268         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
     1269         !  
     1270         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs ) 
     1271         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw ) 
     1272         ! 
    7381273         !               ! ------------------------------------- ! 
    739       CASE( 'F' )        ! interpolation from U-point to F-point ! 
     1274      CASE( 'F' )        ! interpolation from T-point to F-point ! 
    7401275         !               ! ------------------------------------- ! 
    7411276         ! horizontal surface weighted interpolation 
    742          DO jk = 1, jpk 
     1277         IF (nmet==0) THEN 
     1278            DO jk=1,jpk 
     1279               DO jj = 1, jpjm1 
     1280                  DO ji = 1, fs_jpim1   ! vector opt. 
     1281                     pe3_out(ji,jj,jk) = 0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 
     1282                           &      * (  e12t(ji  ,jj  ) * ( pe3_in(ji  ,jj  ,jk)-e3t_0(ji  ,jj  ,jk) )   &  
     1283                           &         + e12t(ji  ,jj+1) * ( pe3_in(ji  ,jj+1,jk)-e3t_0(ji  ,jj+1,jk) )   & 
     1284                           &         + e12t(ji+1,jj  ) * ( pe3_in(ji+1,jj  ,jk)-e3t_0(ji+1,jj  ,jk) )   &  
     1285                           &         + e12t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ))                                                       
     1286                  END DO 
     1287               END DO 
     1288            END DO  
     1289         ENDIF 
     1290         ! 
     1291         IF ( (nmet==1).OR.(nmet==2) ) THEN 
    7431292            DO jj = 1, jpjm1 
    7441293               DO ji = 1, fs_jpim1   ! vector opt. 
    745                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               & 
    746                      &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    747                      &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    748                END DO 
    749             END DO 
    750          END DO 
     1294                  ! bottom correction: 
     1295                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1))  
     1296                  pe3_out(ji,jj,jkbot) = -0.25_wp * umask(ji,jj,jkbot) * umask(ji,jj+1,jkbot) * r1_e12f(ji,jj) & 
     1297                        &                         * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jkbot)   & 
     1298                        &                            + e12t(ji+1,jj  ) * zw(ji+1,jj  ,jkbot)   & 
     1299                        &                            + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jkbot)   & 
     1300                        &                            + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jkbot)   ) 
     1301                  ! 
     1302                  ! If there is a step, taper bottom interface: 
     1303                  IF (hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ) THEN 
     1304                     IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN 
     1305                        IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN 
     1306!                           zmin = ztap * pe3_in(ji+1,jj+1,jkbot) 
     1307                           zmin = ztap * (-zw(ji+1,jj+1,jkbot)+e3t_0(ji+1,jj+1,jkbot)) 
     1308                        ELSE 
     1309!                           zmin = ztap * pe3_in(ji  ,jj+1,jkbot) 
     1310                           zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot)) 
     1311                        ENDIF 
     1312                     ELSE 
     1313                        IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN 
     1314!                           zmin = ztap * pe3_in(ji+1,jj  ,jkbot) 
     1315                           zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot)) 
     1316                        ELSE 
     1317!                           zmin = ztap * pe3_in(ji  ,jj  ,jkbot) 
     1318                           zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot)) 
     1319                        ENDIF 
     1320                     ENDIF 
     1321                     zmin = -e3f_0(ji,jj,jkbot) + zmin 
     1322                     pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin) 
     1323                  ENDIF 
     1324                  ! 
     1325                  zdo =  -pe3_out(ji,jj,jkbot) 
     1326                  DO jk=jkbot-1,1,-1 
     1327                     zup =  0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 
     1328                           &        * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  &  
     1329                           &           + e12t(ji+1,jj  ) * zw(ji+1,jj  ,jk)  &  
     1330                           &           + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  &  
     1331                           &           + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  ) 
     1332                     pe3_out(ji,jj,jk) = zdo - zup  
     1333                     ! 
     1334                     zdo =  zdo - pe3_out(ji,jj,jk)                                     
     1335                  END DO 
     1336                  ! 
     1337               END DO 
     1338            END DO  
     1339         ENDIF 
     1340         ! 
     1341         IF (nmet==2) THEN           ! Spread sea level anomaly 
     1342            ! 
     1343            DO jj = 1, jpjm1 
     1344               DO ji = 1, fs_jpim1   ! vector opt. 
     1345                  DO jk=1,jpk 
     1346                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                         & 
     1347                           &               + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) )                     &  
     1348                           &               / ( hf_0(ji,jj) + 1._wp - umask_i(ji,jj)*umask_i(ji,jj+1) )   & 
     1349                           &               * 0.25_wp * umask(ji,jj,jk)*umask(ji,jj+1,jk)*r1_e12f(ji,jj)  & 
     1350                           &               * ( e12t(ji  ,jj)*zs(ji  ,jj) + e12t(ji  ,jj+1)*zs(ji  ,jj+1) & 
     1351                           &                  +e12t(ji+1,jj)*zs(ji+1,jj) + e12t(ji+1,jj+1)*zs(ji+1,jj+1) )                
     1352                  END DO 
     1353               END DO 
     1354            END DO 
     1355         END IF 
    7511356         ! 
    7521357         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 
     
    7541359         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 
    7551360         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     1361         !  
     1362         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs ) 
     1363         IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw ) 
     1364         ! 
    7561365         !               ! ------------------------------------- ! 
    7571366      CASE( 'W' )        ! interpolation from T-point to W-point ! 
     
    8081417      !! * Local declarations 
    8091418      INTEGER ::   jk 
    810       INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
     1419      INTEGER, DIMENSION(7) ::   id            ! local integers 
     1420      REAL(wp), POINTER, DIMENSION(:,:) :: zhdiv 
    8111421      !!---------------------------------------------------------------------- 
    8121422      ! 
     
    8181428            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    ) 
    8191429            ! 
    820             id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 
    821             id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 
    822             id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
    823             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    824             id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     1430            id(1) = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 
     1431            id(2) = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 
     1432            id(3) = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
     1433            id(4) = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
     1434            id(5) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. ) 
     1435            id(6) = iom_varid( numror, 'un_lf'   , ldstop = .FALSE. ) 
     1436            id(7) = iom_varid( numror, 'vn_lf'   , ldstop = .FALSE. ) 
     1437 
    8251438            !                             ! --------- ! 
    8261439            !                             ! all cases ! 
    8271440            !                             ! --------- ! 
    828             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
     1441            IF( MIN( id(1), id(2) ) > 0 ) THEN       ! all required arrays exist 
    8291442               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    8301443               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 
     
    8381451                  fse3t_b(:,:,:) = fse3t_n(:,:,:) 
    8391452               ENDIF 
    840             ELSE IF( id1 > 0 ) THEN 
     1453            ELSE IF( id(1) > 0 ) THEN 
    8411454               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files' 
    8421455               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.' 
    8431456               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    8441457               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    845                fse3t_n(:,:,:) = fse3t_b(:,:,:) 
     1458               fse3t_b(:,:,:) = fse3t_n(:,:,:) 
    8461459               neuler = 0 
    847             ELSE IF( id2 > 0 ) THEN 
     1460            ELSE IF( id(2) > 0 ) THEN 
    8481461               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files' 
    8491462               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.' 
     
    8671480            IF( ln_vvl_zstar ) THEN       ! z_star case ! 
    8681481               !                          ! ----------- ! 
    869                IF( MIN( id3, id4 ) > 0 ) THEN 
     1482               IF( MIN( id(3), id(4) ) > 0 ) THEN 
    8701483                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 
    8711484               ENDIF 
     
    8731486            ELSE                          ! z_tilde and layer cases ! 
    8741487               !                          ! ----------------------- ! 
    875                IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist 
     1488               IF( MIN( id(3), id(4) ) > 0 ) THEN  ! all required arrays exist 
    8761489                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 
    8771490                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) ) 
     
    8791492                  tilde_e3t_b(:,:,:) = 0.0_wp 
    8801493                  tilde_e3t_n(:,:,:) = 0.0_wp 
    881                ENDIF 
    882                !                          ! ------------ ! 
     1494                  ! 
     1495                  neuler = 0 
     1496               ENDIF                      ! ------------ ! 
    8831497               IF( ln_vvl_ztilde ) THEN   ! z_tilde case ! 
    8841498                  !                       ! ------------ ! 
    885                   IF( id5 > 0 ) THEN  ! required array exists 
    886                      CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) ) 
    887                   ELSE                ! array is missing 
    888                      hdiv_lf(:,:,:) = 0.0_wp 
     1499                  IF( MINVAL(id(5:7) ) > 0 ) THEN  ! all required arrays exist 
     1500                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) ) 
     1501                     CALL iom_get( numror, jpdom_autoglo, 'un_lf'   ,    un_lf(:,:,:) ) 
     1502                     CALL iom_get( numror, jpdom_autoglo, 'vn_lf'   ,    vn_lf(:,:,:) ) 
     1503                  ELSE                            ! one at least array is missing 
     1504                     hdivn_lf(:,:,:) = 0.0_wp 
     1505                     un_lf(:,:,:)    = 0.0_wp 
     1506                     vn_lf(:,:,:)    = 0.0_wp 
     1507                     neuler = 0 
    8891508                  ENDIF 
    8901509               ENDIF 
     1510               ! 
    8911511            ENDIF 
    8921512            ! 
     
    8981518               tilde_e3t_b(:,:,:) = 0.0_wp 
    8991519               tilde_e3t_n(:,:,:) = 0.0_wp 
    900                IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp 
    9011520            END IF 
     1521            IF( ln_vvl_ztilde ) THEN 
     1522               hdivn_lf(:,:,:) = 0.0_wp 
     1523               un_lf(:,:,:)    = 0.0_wp 
     1524               vn_lf(:,:,:)    = 0.0_wp 
     1525            ENDIF 
    9021526         ENDIF 
    9031527 
     
    9191543         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    9201544            !                                        ! ------------ ! 
    921             CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 
     1545            CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:) ) 
     1546            CALL iom_rstput( kt, nitrst, numrow, 'un_lf'   , un_lf(:,:,:)    ) 
     1547            CALL iom_rstput( kt, nitrst, numrow, 'vn_lf'   , vn_lf(:,:,:)    ) 
    9221548         ENDIF 
    9231549 
     
    9261552 
    9271553   END SUBROUTINE dom_vvl_rst 
    928  
    9291554 
    9301555   SUBROUTINE dom_vvl_ctl 
     
    9381563      INTEGER ::   ios 
    9391564 
    940       NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 
    941                       & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , & 
    942                       & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe 
     1565      NAMELIST/nam_vvl/ ln_vvl_zstar               , ln_vvl_ztilde                       , & 
     1566                      & ln_vvl_layer               , ln_vvl_ztilde_as_zstar              , & 
     1567                      & ln_vvl_zstar_at_eqtor      , ln_vvl_zstar_on_shelf               , & 
     1568                      & ln_vvl_adv_cn2             , ln_vvl_adv_fct                      , & 
     1569                      & ln_vvl_lap                 , ln_vvl_blp                          , & 
     1570                      & rn_ahe3_lap                , rn_ahe3_blp                         , & 
     1571                      & rn_rst_e3t                 , rn_lf_cutoff                        , & 
     1572                      & ln_vvl_regrid              , rn_zdef_max                         , & 
     1573                      & ln_vvl_dbg   ! not yet implemented: ln_vvl_kepe 
    9431574      !!----------------------------------------------------------------------  
    9441575 
     
    9611592         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer 
    9621593         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar 
    963          WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor 
    9641594         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation' 
    9651595         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe 
    966          WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient' 
    967          WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3 
     1596         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor 
     1597         WRITE(numout,*) '      ztilde on shelves          ln_vvl_zstar_on_shelf  = ', ln_vvl_zstar_on_shelf 
     1598         WRITE(numout,*) '           Namelist nam_vvl : thickness advection scheme' 
     1599         WRITE(numout,*) '              2nd order                  ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2 
     1600         WRITE(numout,*) '              2nd order FCT              ln_vvl_adv_fct = ', ln_vvl_adv_fct                     
     1601         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion scheme' 
     1602         WRITE(numout,*) '              Laplacian                  ln_vvl_lap     = ', ln_vvl_lap 
     1603         WRITE(numout,*) '              Bilaplacian                ln_vvl_blp     = ', ln_vvl_blp 
     1604         WRITE(numout,*) '              Laplacian   coefficient    rn_ahe3_lap    = ', rn_ahe3_lap 
     1605         WRITE(numout,*) '              Bilaplacian coefficient    rn_ahe3_blp    = ', rn_ahe3_blp 
     1606         WRITE(numout,*) '           Namelist nam_vvl : layers regriding' 
     1607         WRITE(numout,*) '              ln_vvl_regrid                              = ', ln_vvl_regrid 
    9681608         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change' 
    9691609         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max 
     
    9921632 
    9931633      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    994       IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
     1634 
     1635      IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN  
     1636         ioptio = 0                      ! Choose one advection scheme at most 
     1637         IF( ln_vvl_adv_cn2         )        ioptio = ioptio + 1 
     1638         IF( ln_vvl_adv_fct         )        ioptio = ioptio + 1 
     1639         IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' ) 
     1640      ENDIF 
    9951641 
    9961642      IF(lwp) THEN                   ! Print the choice 
     
    14012047   END SUBROUTINE dom_vvl_orca_fix 
    14022048 
     2049   SUBROUTINE dom_vvl_zdf( kt, p2dt )  
     2050      !!---------------------------------------------------------------------- 
     2051      !!                  ***  ROUTINE dom_vvl_zdf  *** 
     2052      !! 
     2053      !! ** Purpose :   Do vertical thicknesses anomaly diffusion 
     2054      !! 
     2055      !! ** Method  :   
     2056      !! 
     2057      !! ** Action  :  
     2058      !!--------------------------------------------------------------------- 
     2059      USE oce     , ONLY:   zwd => ua       , zws => va         ! (ua,va) used as 3D workspace 
     2060      ! 
     2061      INTEGER                         , INTENT(in) ::   kt     ! ocean time-step index 
     2062      REAL(wp)                        , INTENT(in) ::   p2dt   ! time step 
     2063      ! 
     2064      INTEGER  :: ji, jj, jk ! dummy loop indices 
     2065      REAL(wp) :: zr_tscale 
     2066      REAL(wp) :: za1, za2, za3, za4, zdiff 
     2067      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt 
     2068      !!--------------------------------------------------------------------- 
     2069      ! 
     2070      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_zdf') 
     2071      ! 
     2072      zr_tscale = 0.5 
     2073      ! 
     2074      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt)  
     2075      ! 
     2076      ! 
     2077 
     2078      zwt(:,:,1:jpk) = 1.e-10 
     2079 
     2080      zwt(:,:,1) = 0.0_wp 
     2081!      DO jk = 2, jpkm1 
     2082!         DO jj = 1, jpj 
     2083!            DO ji = 1, jpi 
     2084!               zwt(ji,jj,jk) = zwt(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     2085!            END DO 
     2086!         END DO 
     2087!      END DO 
     2088      ! 
     2089      ! 
     2090      ! Set diffusivity (homogeneous to an inverse time scale) 
     2091      ! 
     2092      DO jk = 2, jpkm1 
     2093         DO jj = 2, jpjm1 
     2094            DO ji = fs_2, fs_jpim1 
     2095! Taper a little bit: 
     2096                za1 = tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1) 
     2097                za2 = tilde_e3t_n(ji,jj,jk  )+e3t_0(ji,jj,jk  ) 
     2098                za4 = 0.5_wp * (e3t_0(ji,jj,jk-1) + e3t_0(ji,jj,jk  )) 
     2099                za3 = 0.5_wp * (za1 + za2)  
     2100                zdiff = ABS(za3-za4)/za4 
     2101!                IF (zdiff>=0.8) THEN 
     2102!                   zwt(ji,jj,jk) =  zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk) 
     2103                   zwt(ji,jj,jk) =  dsm(ji,jj)/ht_0(ji,jj)*(1._wp-tanh((mbkt(ji,jj)+1-jk)*0.2))*tmask(ji,jj,jk) 
     2104 
     2105!                ELSE 
     2106!                   zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk) 
     2107!                ENDIF 
     2108            END DO 
     2109         END DO 
     2110      END DO 
     2111      !  
     2112      ! 
     2113      ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
     2114      DO jk = 1, jpkm1 
     2115         DO jj = 2, jpjm1 
     2116            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2117                zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / fse3w(ji,jj,jk  ) 
     2118                zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / fse3w(ji,jj,jk+1) 
     2119                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     2120            END DO 
     2121         END DO 
     2122      END DO 
     2123      ! 
     2124      !! Matrix inversion from the first level 
     2125      !!---------------------------------------------------------------------- 
     2126      DO jj = 2, jpjm1 
     2127         DO ji = fs_2, fs_jpim1 
     2128            zwt(ji,jj,1) = zwd(ji,jj,1) 
     2129         END DO 
     2130      END DO 
     2131      DO jk = 2, jpkm1 
     2132         DO jj = 2, jpjm1 
     2133            DO ji = fs_2, fs_jpim1 
     2134               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
     2135            END DO 
     2136         END DO 
     2137      END DO 
     2138      !          
     2139      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
     2140      DO jk = 2, jpkm1 
     2141         DO jj = 2, jpjm1 
     2142            DO ji = fs_2, fs_jpim1 
     2143               tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * tilde_e3t_a(ji,jj,jk-1) 
     2144            END DO 
     2145         END DO 
     2146      END DO 
     2147      ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
     2148      DO jj = 2, jpjm1 
     2149         DO ji = fs_2, fs_jpim1 
     2150            tilde_e3t_a(ji,jj,jpkm1) = tilde_e3t_a(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     2151         END DO 
     2152      END DO 
     2153      DO jk = jpk-2, 1, -1 
     2154         DO jj = 2, jpjm1 
     2155            DO ji = fs_2, fs_jpim1 
     2156               tilde_e3t_a(ji,jj,jk) = ( tilde_e3t_a(ji,jj,jk) - zws(ji,jj,jk) * tilde_e3t_a(ji,jj,jk+1) )   & 
     2157                  &            / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
     2158            END DO 
     2159         END DO 
     2160      END DO 
     2161      ! 
     2162      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt)  
     2163      ! 
     2164      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_zdf') 
     2165      ! 
     2166   END SUBROUTINE dom_vvl_zdf 
     2167 
     2168   SUBROUTINE dom_vvl_zdf2( kt, p2dt )  
     2169      !!---------------------------------------------------------------------- 
     2170      !!                  ***  ROUTINE dom_vvl_zdf  *** 
     2171      !! 
     2172      !! ** Purpose :   Do vertical interface diffusion 
     2173      !! 
     2174      !! ** Method  :   
     2175      !! 
     2176      !! ** Action  :  
     2177      !!--------------------------------------------------------------------- 
     2178      USE oce     , ONLY:   zwd => ua       , zws => va         ! (ua,va) used as 3D workspace 
     2179      ! 
     2180      INTEGER                         , INTENT(in) ::   kt     ! ocean time-step index 
     2181      REAL(wp)                        , INTENT(in) ::   p2dt   ! time step 
     2182      ! 
     2183      INTEGER  :: ji, jj, jk, kbot, kbotm1 ! dummy loop indices 
     2184      REAL(wp) :: zr_tscale 
     2185      REAL(wp) :: za1, za2, za3, za4, zdiff 
     2186      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt, zw 
     2187      !!--------------------------------------------------------------------- 
     2188      ! 
     2189      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_zdf') 
     2190      ! 
     2191      zr_tscale = 0.5 
     2192      ! 
     2193      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zw)  
     2194      ! 
     2195      ! Compute internal interfaces depths: 
     2196      zw(:,:,1) = 0._wp 
     2197      DO jk = 2, jpk 
     2198         DO jj = 2, jpjm1 
     2199            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2200               zw(ji,jj,jk) = zw(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))  * tmask(ji,jj,jk-1) 
     2201            END DO 
     2202         END DO 
     2203      END DO 
     2204      ! 
     2205      ! Set diffusivities at interfaces 
     2206      zwt(:,:,:) = 0.00000001_wp * tmask(:,:,:)       
     2207      zwt(:,:,1) = 0._wp 
     2208      !  
     2209      ! 
     2210      ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
     2211      DO jk = 2, jpkm1 
     2212         DO jj = 2, jpjm1 
     2213            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2214                zwi(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk-1) + zwt(ji,jj,jk  )) &  
     2215                              &                   / fse3w(ji,jj,jk-1) / fse3t(ji,jj,jk-1) & 
     2216                              &                   * ht(ji,jj)  
     2217                zws(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk  ) + zwt(ji,jj,jk+1)) &  
     2218                              &                   / fse3w(ji,jj,jk  ) / fse3t(ji,jj,jk  ) & 
     2219                              &                   * ht(ji,jj)  
     2220 
     2221                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     2222            END DO 
     2223         END DO 
     2224      END DO 
     2225      ! Boundary conditions (Neumann) 
     2226      DO jj = 2, jpjm1 
     2227         DO ji = fs_2, fs_jpim1 
     2228            zwi(ji,jj,1) = 0._wp 
     2229            zws(ji,jj,1) = 0._wp 
     2230            zwd(ji,jj,1) = 1._wp 
     2231            zw (ji,jj,1) = 0._wp 
     2232            ! 
     2233            zwd(ji,jj,2) = zwd(ji,jj,2) +  zwi(ji,jj,2) 
     2234            zwi(ji,jj,2) = 0._wp 
     2235!            zwi(ji,jj,2) = 0._wp 
     2236!            zws(ji,jj,2) = 0._wp 
     2237!            zwd(ji,jj,2) = 1._wp 
     2238         END DO 
     2239      END DO 
     2240      DO jj = 2, jpjm1 
     2241         DO ji = fs_2, fs_jpim1 
     2242            kbot   = mbkt(ji,jj) + 1 
     2243            kbotm1 = mbkt(ji,jj) 
     2244            zwi(ji,jj,kbot  ) = 0._wp 
     2245            zws(ji,jj,kbot  ) = 0._wp 
     2246            zwd(ji,jj,kbot  ) = 1._wp 
     2247            ! 
     2248            zwd(ji,jj,kbotm1) = zwd(ji,jj,kbotm1) +  zws(ji,jj,kbotm1) 
     2249            zws(ji,jj,kbotm1) = 0._wp 
     2250!            zwi(ji,jj,kbotm1) = 0._wp 
     2251!            zws(ji,jj,kbotm1) = 0._wp 
     2252!            zwd(ji,jj,kbotm1) = 1._wp 
     2253         END DO 
     2254      END DO 
     2255      ! 
     2256      DO jk = 2, jpkm1 
     2257         DO jj = 2, jpjm1 
     2258            DO ji = fs_2, fs_jpim1 
     2259               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
     2260            END DO 
     2261         END DO 
     2262      END DO 
     2263      ! 
     2264      DO jk = 2, jpk 
     2265         DO jj = 2, jpjm1 
     2266            DO ji = fs_2, fs_jpim1    ! vector opt. 
     2267               zwi(ji,jj,jk) = zw(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) *zwi(ji,jj,jk-1) 
     2268            END DO 
     2269         END DO 
     2270      END DO 
     2271      ! 
     2272      DO jk = jpk-1, 2, -1 
     2273         DO jj = 2, jpjm1 
     2274            DO ji = fs_2, fs_jpim1    ! vector opt. 
     2275               zw(ji,jj,jk) = ( zwi(ji,jj,jk) - zws(ji,jj,jk) * zw(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     2276            END DO 
     2277         END DO 
     2278      END DO 
     2279      ! 
     2280      ! Revert to thicknesses anomalies: 
     2281      DO jk = 1, jpkm1 
     2282         DO jj = 2, jpjm1 
     2283            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2284               tilde_e3t_a(ji,jj,jk) = (zw(ji,jj,jk+1)-zw(ji,jj,jk)-e3t_0(ji,jj,jk))* tmask(ji,jj,jk) 
     2285            END DO 
     2286         END DO 
     2287      END DO 
     2288      ! 
     2289      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zw)  
     2290      ! 
     2291      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_zdf') 
     2292      ! 
     2293   END SUBROUTINE dom_vvl_zdf2 
     2294 
     2295   SUBROUTINE dom_vvl_regrid( kt ) 
     2296      !!---------------------------------------------------------------------- 
     2297      !!                ***  ROUTINE dom_vvl_regrid  *** 
     2298      !!                    
     2299      !! ** Purpose :  Ensure "well-behaved" vertical grid 
     2300      !! 
     2301      !! ** Method  :  More or less adapted from references below. 
     2302      !! 
     2303      !! ** Action  :  Ensure that thickness are above a given value, spaced enough 
     2304      !!               and revert to Eulerian coordinates near the bottom.       
     2305      !! 
     2306      !! References :  Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction 
     2307      !!               with a Model Combining Terrain-following and Isentropic 
     2308      !!               coordinates. Part I: Model Description. Monthly Weather Rev., 
     2309      !!               121, 1770-1785. 
     2310      !!               Toy, M., 2011: Incorporating Condensational Heating into a 
     2311      !!               Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic- 
     2312      !!               Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954. 
     2313      !!---------------------------------------------------------------------- 
     2314      !! * Arguments 
     2315      INTEGER, INTENT( in )               :: kt         ! time step 
     2316 
     2317      !! * Local declarations 
     2318      INTEGER                             :: ji, jj, jk ! dummy loop indices 
     2319      LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond 
     2320      LOGICAL                             :: ll_e3tsurf_const, ll_zdiff_cond, ll_blpdiff_cond 
     2321      INTEGER                             :: jkbot 
     2322      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd, dzmin_vvl 
     2323      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj 
     2324      REAL(wp)                            :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn 
     2325      REAL(wp)                            :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot 
     2326      REAL(wp)                            :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim 
     2327      REAL(wp), POINTER, DIMENSION(:,:)   :: zdw 
     2328      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zwdw_b 
     2329      !!---------------------------------------------------------------------- 
     2330 
     2331      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_regrid') 
     2332      ! 
     2333      CALL wrk_alloc( jpi, jpj, jpk, zwdw)  
     2334      ! 
     2335      ! Some user defined parameters below: 
     2336      ll_chk_bot2top   = .TRUE. 
     2337      ll_chk_top2bot   = .TRUE. 
     2338      ll_e3tsurf_const = .FALSE. 
     2339      dzmin_vvl  = 1._wp   ! Absolute minimum depth (in meters) 
     2340      zfrch_stp  = 0.5_wp   ! Maximum fractionnal thickness change in one time step (<= 1.) 
     2341      zfrch_rel  = 0.0_wp   ! Maximum relative thickness change in the vertical (<= 1.) 
     2342      zfrac_bot  = 0.01_wp  ! Fraction of bottom level allowed to change 
     2343      zscal_bot  = 2._wp    ! Depth lengthscale 
     2344      ll_zdiff_cond    = .FALSE.  ! Conditionnal vertical diffusion of interfaces 
     2345         zvdiff  = 0.1_wp   ! m/s 
     2346         zvlim   = 0.4_wp   ! max d2h/dh 
     2347      ll_lapdiff_cond  = .FALSE.  ! Conditionnal Laplacian diffusion of interfaces 
     2348         zhdiff  = 0.1_wp   ! ad. 
     2349         zhlim   = 0.01_wp  ! max lap(z)/e1 
     2350      ll_blpdiff_cond  = .FALSE.  ! Conditionnal Bilaplacian diffusion of interfaces 
     2351         zhdiff2 = 0.8_wp   ! ad. 
     2352         zhlim2  = 5.e-11_wp  ! max bilap(z)/e1**3 
     2353      ! ---------------------------------------------------------------------------------------  
     2354      ! 
     2355      ! Set arrays determining maximum vertical displacement at the bottom: 
     2356      !-------------------------------------------------------------------- 
     2357      IF ( kt==nit000 ) THEN 
     2358         DO jj = 2, jpjm1 
     2359            DO ji = 2, jpim1 
     2360               jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1)) 
     2361               i_int_bot(ji,jj) = jk 
     2362            END DO 
     2363         END DO 
     2364         dsm(:,:) = REAL( i_int_bot(:,:), wp )   ;   CALL lbc_lnk(dsm(:,:),'T',1.) 
     2365         i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 ) 
     2366 
     2367         CALL wrk_alloc( jpi, jpj, zdw ) 
     2368         DO jj = 2, jpjm1 
     2369            DO ji = 2, jpim1 
     2370               zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji  ,jj,1), & 
     2371                          &     ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), & 
     2372                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj  ,1), & 
     2373                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1)  ) 
     2374               zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall )  
     2375            END DO 
     2376         END DO 
     2377         CALL lbc_lnk( zdw(:,:), 'T', 1. ) 
     2378 
     2379         DO jj = 2, jpjm1 
     2380            DO ji = 2, jpim1 
     2381               dsm(ji,jj) = 1._wp/16._wp * (            zdw(ji-1,jj-1) + zdw(ji+1,jj-1)      & 
     2382                  &                           +         zdw(ji-1,jj+1) + zdw(ji+1,jj+1)      & 
     2383                  &                           + 2._wp*( zdw(ji  ,jj-1) + zdw(ji-1,jj  )      & 
     2384                  &                           +         zdw(ji+1,jj  ) + zdw(ji  ,jj+1) )    & 
     2385                  &                           + 4._wp*  zdw(ji  ,jj  )                       ) 
     2386            END DO 
     2387         END DO 
     2388 
     2389         CALL lbc_lnk( dsm(:,:), 'T', 1. ) 
     2390         CALL wrk_dealloc( jpi, jpj, zdw)          
     2391      
     2392         IF (ln_zps) THEN 
     2393            DO jj = 1, jpj 
     2394               DO ji = 1, jpi 
     2395                  jk = i_int_bot(ji,jj) 
     2396                  hsm(ji,jj) = zfrac_bot * e3w_1d(jk) 
     2397               END DO 
     2398            END DO 
     2399         ELSE 
     2400            DO jj = 1, jpj 
     2401               DO ji = 1, jpi 
     2402                  jk = i_int_bot(ji,jj) 
     2403                  hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk) 
     2404               END DO 
     2405            END DO 
     2406         ENDIF 
     2407      END IF 
     2408 
     2409      ! Provisionnal interface depths: 
     2410      !------------------------------- 
     2411      zwdw(:,:,1) = 0.e0 
     2412      DO jj = 1, jpj 
     2413         DO ji = 1, jpi 
     2414            DO jk = 2, jpk 
     2415               zwdw(ji,jj,jk) =  zwdw(ji,jj,jk-1) + &  
     2416                              & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     2417            END DO 
     2418         END DO 
     2419      END DO 
     2420      ! 
     2421      ! Conditionnal horizontal Laplacian diffusion: 
     2422      !--------------------------------------------- 
     2423      IF ( ll_lapdiff_cond ) THEN 
     2424         CALL wrk_alloc( jpi, jpj, jpk, zwdw_b) 
     2425         ! 
     2426         zwdw_b(:,:,1) = 0._wp 
     2427         DO jj = 1, jpj 
     2428            DO ji = 1, jpi 
     2429               DO jk=2,jpk 
     2430                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + &  
     2431                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     2432               END DO 
     2433            END DO 
     2434         END DO 
     2435         ! 
     2436         DO jk = 2, jpkm1 
     2437            DO jj = 1, jpjm1 
     2438               DO ji = 1, fs_jpim1   ! vector opt.                   
     2439                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     2440                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) ) 
     2441                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     2442                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) ) 
     2443               END DO 
     2444            END DO 
     2445         END DO 
     2446             
     2447         DO jk = 2, jpkm1 
     2448            DO jj = 2, jpjm1 
     2449               DO ji = fs_2, fs_jpim1   ! vector opt. 
     2450                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    & 
     2451                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj) 
     2452                  zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e12t(ji,jj)), 0._wp) 
     2453                  ztmp = SIGN(zh2, ztmp1) 
     2454                  zeu2 = zhdiff * e12t(ji,jj)*e12t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj)) 
     2455                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk) 
     2456               END DO 
     2457            END DO 
     2458         END DO          
     2459         ! 
     2460         CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b) 
     2461         ! 
     2462      ENDIF 
     2463 
     2464      ! Conditionnal horizontal Bilaplacian diffusion: 
     2465      !----------------------------------------------- 
     2466      IF ( ll_blpdiff_cond ) THEN 
     2467         CALL wrk_alloc( jpi, jpj, jpk, zwdw_b) 
     2468         ! 
     2469         zwdw_b(:,:,1) = 0._wp 
     2470         DO jj = 1, jpj 
     2471            DO ji = 1, jpi 
     2472               DO jk=2,jpk 
     2473                  zwdw_b(ji,jj,jk) =  zwdw_b(ji,jj,jk-1) + &  
     2474                                   & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 
     2475               END DO 
     2476            END DO 
     2477         END DO 
     2478         ! 
     2479         DO jk = 2, jpkm1 
     2480            DO jj = 1, jpjm1 
     2481               DO ji = 1, fs_jpim1   ! vector opt.                   
     2482                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     2483                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) ) 
     2484                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     2485                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) ) 
     2486               END DO 
     2487            END DO 
     2488         END DO 
     2489             
     2490         DO jk = 2, jpkm1 
     2491            DO jj = 2, jpjm1 
     2492               DO ji = fs_2, fs_jpim1   ! vector opt. 
     2493                  zwdw_b(ji,jj,jk) = -( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    & 
     2494                                &  +    (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj) 
     2495               END DO 
     2496            END DO 
     2497         END DO    
     2498         ! 
     2499         CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )       
     2500         ! 
     2501         DO jk = 2, jpkm1 
     2502            DO jj = 1, jpjm1 
     2503               DO ji = 1, fs_jpim1   ! vector opt.                   
     2504                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     2505                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj  ,jk) ) 
     2506                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     2507                                  &  * ( zwdw_b(ji,jj,jk) - zwdw_b(ji  ,jj+1,jk) ) 
     2508               END DO 
     2509            END DO 
     2510         END DO 
     2511         !    
     2512         DO jk = 2, jpkm1 
     2513            DO jj = 2, jpjm1 
     2514               DO ji = fs_2, fs_jpim1   ! vector opt. 
     2515                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    & 
     2516                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj) 
     2517                  zh2 = MAX(abs(ztmp1)-zhlim2, 0._wp) 
     2518                  ztmp = SIGN(zh2, ztmp1) 
     2519                  zeu2 = zhdiff2 * e12t(ji,jj)*e12t(ji,jj) / 16._wp 
     2520                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp *  tmask(ji,jj,jk) 
     2521               END DO 
     2522            END DO 
     2523         END DO    
     2524         ! 
     2525         CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b) 
     2526         ! 
     2527      ENDIF 
     2528 
     2529      ! Conditionnal vertical diffusion: 
     2530      !--------------------------------- 
     2531      IF ( ll_zdiff_cond ) THEN 
     2532         DO jk = 2, jpkm1 
     2533            DO jj = 2, jpjm1 
     2534               DO ji = fs_2, fs_jpim1   ! vector opt.     
     2535                  ztmp  = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) &  
     2536                            -(tilde_e3t_b(ji,jj,jk  )+e3t_0(ji,jj,jk  ))*tmask(ji,jj,jk  ) )  
     2537                  ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) &  
     2538                        &           +tilde_e3t_b(ji,jj,jk  ) + e3t_0(ji,jj,jk  ) ) 
     2539                  zh2   = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp) 
     2540                  ztmp  = SIGN(zh2, ztmp) 
     2541                  IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0 
     2542                  zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk) 
     2543               END DO 
     2544            END DO 
     2545         END DO 
     2546      ENDIF 
     2547      ! 
     2548      ! Check grid from the bottom to the surface 
     2549      !------------------------------------------ 
     2550      IF ( ll_chk_bot2top ) THEN 
     2551         DO jj = 2, jpjm1 
     2552            DO ji = 2, jpim1 
     2553               jkbot = mbkt(ji,jj)                    
     2554               DO jk = jkbot,2,-1 
     2555                  ! 
     2556                  zh_0   = e3t_0(ji,jj,jk) 
     2557                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)) 
     2558                  zh_old = tilde_e3t_a(ji,jj,jk  ) + zh_0 
     2559!                  zh_dwn = tilde_e3t_a(ji,jj,jk+1) + e3t_0(ji,jj,jk+1) 
     2560                  zh_min = MIN(zh_0/3._wp, dzmin_vvl) 
     2561                  !  
     2562                  ! Set maximum and minimum vertical excursions    
     2563                  ztmph = hsm(ji,jj) 
     2564                  ztmpd = dsm(ji,jj) 
     2565                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)                
     2566                  zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 ) 
     2567                  zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff) 
     2568                  zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 ) 
     2569                  zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff) 
     2570                  ! 
     2571                  ! New layer thickness: 
     2572                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
     2573                  ! 
     2574                  ! Ensure minimum layer thickness:                   
     2575!                  zh_new = MIN(zh_old, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) ) 
     2576                  zh_new = cush(zh_new, zh_min) 
     2577                  ! 
     2578                  ! Final flux: 
     2579                  zdiff = (zh_new - zh_old) * tmask(ji,jj,jk) 
     2580                  ! 
     2581                  ! Limit flux:  
     2582                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 
     2583                  zdiff = SIGN(ztmp, zh_new - zh_old) 
     2584                  zh_new = zdiff + zh_old 
     2585                  ! 
     2586                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
     2587                  zwdw(ji,jj,jk)          = zwdw(ji,jj,jk+1) - zh_new 
     2588                  tilde_e3t_a(ji,jj,jk-1) = (-zdiff + tilde_e3t_a(ji,jj,jk-1) ) * tmask(ji,jj,jk-1)             
     2589               END DO     
     2590            END DO  
     2591         END DO 
     2592      END IF     
     2593      ! 
     2594      ! Check grid from the surface to the bottom  
     2595      !------------------------------------------  
     2596      IF ( ll_chk_top2bot ) THEN       
     2597         DO jj = 2, jpjm1 
     2598            DO ji = 2, jpim1 
     2599               jkbot = mbkt(ji,jj)    
     2600               DO jk = 1, jkbot-1 
     2601                  ! 
     2602                  zh_0   = e3t_0(ji,jj,jk) 
     2603                  zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)) 
     2604                  zh_old = tilde_e3t_a(ji,jj,jk  ) + zh_0    
     2605!                  zh_up  = tilde_e3t_a(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) 
     2606                  zh_min = MIN(zh_0/3._wp, dzmin_vvl) 
     2607                  IF ((jk<=5).AND.ll_e3tsurf_const) zh_min = MAX(e3t_0(ji,jj,1)/3._wp, dzmin_vvl) 
     2608                  ! 
     2609                  ! Ensure minimum layer thickness: 
     2610!                  zh_new=MIN(zh_old, zh_up * zfrch_rel / (2._wp-zfrch_rel) ) 
     2611                  zh_new = cush(zh_old, zh_min) 
     2612                  ! 
     2613                  ! Final flux: 
     2614                  zdiff = (zh_new -zh_old) * tmask(ji,jj,jk) 
     2615                  !  
     2616                  ! Limit flux:                  
     2617                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 
     2618                  zdiff = SIGN(ztmp, zdiff) 
     2619                  zh_new = zdiff + zh_old 
     2620                  ! 
     2621                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
     2622                  zwdw(ji,jj,jk+1)        = zwdw(ji,jj,jk) + zh_new 
     2623                  tilde_e3t_a(ji,jj,jk+1) = (-zdiff + tilde_e3t_a(ji,jj,jk+1) ) * tmask(ji,jj,jk+1) 
     2624               END DO 
     2625               !                
     2626            END DO 
     2627         END DO 
     2628      ENDIF 
     2629      ! 
     2630      CALL wrk_dealloc( jpi, jpj, jpk, zwdw )  
     2631      ! 
     2632      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_regrid') 
     2633      ! 
     2634   END SUBROUTINE dom_vvl_regrid 
     2635 
     2636   FUNCTION cush(hin, hmin)  RESULT(hout) 
     2637      !!---------------------------------------------------------------------- 
     2638      !!                 ***  FUNCTION cush  *** 
     2639      !! 
     2640      !! ** Purpose :    
     2641      !! 
     2642      !! ** Method  : 
     2643      !! 
     2644      !!---------------------------------------------------------------------- 
     2645      IMPLICIT NONE 
     2646      REAL(wp), INTENT(in) ::  hin, hmin 
     2647      REAL(wp)             ::  hout, zx, zh_cri 
     2648      !!---------------------------------------------------------------------- 
     2649      zh_cri = 3._wp * hmin 
     2650      ! 
     2651      IF ( hin<=0._wp ) THEN 
     2652         hout = hmin 
     2653      ! 
     2654      ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN 
     2655         zx = hin/zh_cri 
     2656         hout = hmin * (1._wp + zx + zx*zx)       
     2657      ! 
     2658      ELSEIF ( hin>zh_cri ) THEN 
     2659         hout = hin 
     2660      ! 
     2661      ENDIF 
     2662      ! 
     2663   END FUNCTION cush 
     2664 
     2665   FUNCTION cush_max(hin, hmax)  RESULT(hout) 
     2666      !!---------------------------------------------------------------------- 
     2667      !!                 ***  FUNCTION cush  *** 
     2668      !! 
     2669      !! ** Purpose :    
     2670      !! 
     2671      !! ** Method  : 
     2672      !! 
     2673      !!---------------------------------------------------------------------- 
     2674      IMPLICIT NONE 
     2675      REAL(wp), INTENT(in) ::  hin, hmax 
     2676      REAL(wp)             ::  hout, hmin, zx, zh_cri 
     2677      !!---------------------------------------------------------------------- 
     2678      hmin = 0.1_wp * hmax 
     2679      zh_cri = 3._wp * hmin 
     2680      ! 
     2681      IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN 
     2682         zx = (hmax-hin)/zh_cri 
     2683         hout = hmax - hmin * (1._wp + zx + zx*zx) 
     2684         ! 
     2685      ELSEIF ( hin>(hmax-zh_cri) ) THEN 
     2686         hout = hmax - hmin 
     2687         ! 
     2688      ELSE 
     2689         hout = hin 
     2690         ! 
     2691      ENDIF 
     2692      ! 
     2693   END FUNCTION cush_max 
     2694 
     2695   SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin ) 
     2696      !!---------------------------------------------------------------------- 
     2697      !!                  ***  ROUTINE dom_vvl_adv_fct  *** 
     2698      !!  
     2699      !! **  Purpose :  Do thickness advection 
     2700      !! 
     2701      !! **  Method  :  FCT scheme to ensure positivity  
     2702      !! 
     2703      !! **  Action  : - Update pta thickness tendency and diffusive fluxes 
     2704      !!               - this is the total trend, hence it does include sea level motions 
     2705      !!               - Upstream corrections to antidiffusive fluxes ensure 
     2706      !!                 that barotropic transport matches what is contained in input fluxes 
     2707      !!---------------------------------------------------------------------- 
     2708      ! 
     2709      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     2710      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta    ! thickness baroclinic trend  
     2711      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)    ::   uin, vin  ! input velocities 
     2712      ! 
     2713      INTEGER  ::   ji, jj, jk               ! dummy loop indices   
     2714      INTEGER  ::   ikbu, ikbv, ibot 
     2715      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar 
     2716      REAL(wp) ::   zdi, zdj, zmin           !   -      - 
     2717      REAL(wp) ::   zfp_ui, zfp_vj           !   -      - 
     2718      REAL(wp) ::   zfm_ui, zfm_vj           !   -      - 
     2719      REAL(wp) ::   zfp_hi, zfp_hj           !   -      - 
     2720      REAL(wp) ::   zfm_hi, zfm_hj           !   -      - 
     2721      REAL(wp), POINTER, DIMENSION(:,:)   :: zbu, zbv, zhu_b, zhv_b 
     2722      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwi 
     2723      !!---------------------------------------------------------------------- 
     2724      ! 
     2725      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_adv_fct') 
     2726      ! 
     2727      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv) 
     2728      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwi) 
     2729      ! 
     2730      ! 1. Initializations 
     2731      ! ------------------ 
     2732      ! 
     2733      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     2734         z2dtt =  rdt 
     2735      ELSE 
     2736         z2dtt = 2.0_wp * rdt 
     2737      ENDIF 
     2738      ! 
     2739      zwi(:,:,:) = 0.e0 
     2740      zwx(:,:,:) = 0.e0  
     2741      zwy(:,:,:) = 0.e0 
     2742      ! 
     2743      ! 
     2744      ! 2. upstream advection with initial mass fluxes & intermediate update 
     2745      ! -------------------------------------------------------------------- 
     2746      DO jk = 1, jpkm1 
     2747         DO jj = 1, jpjm1 
     2748            DO ji = 1, fs_jpim1   ! vector opt. 
     2749               ! 
     2750               zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 
     2751               zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 
     2752               zfp_hi = fse3t_b(ji  ,jj  ,jk) 
     2753               zfm_hi = fse3t_b(ji+1,jj  ,jk) 
     2754               zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2755               ! 
     2756               zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 
     2757               zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 
     2758               zfp_hj = fse3t_b(ji  ,jj  ,jk) 
     2759               zfm_hj = fse3t_b(ji  ,jj+1,jk) 
     2760               zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2761            END DO 
     2762         END DO 
     2763      END DO 
     2764 
     2765      IF ( .NOT.ln_sco ) THEN 
     2766         ! Correct bottom upstream fluxes 
     2767         ! considering "shelf horizon depths" 
     2768         DO jj = 1, jpjm1 
     2769            DO ji = 1, fs_jpim1   ! vector opt. 
     2770               ! upstream scheme 
     2771               ikbu = mbku(ji,jj) 
     2772               ikbv = mbkv(ji,jj) 
     2773               zfp_ui = uin(ji,jj,ikbu) + ABS( uin(ji,jj,ikbu) ) 
     2774               zfm_ui = uin(ji,jj,ikbu) - ABS( uin(ji,jj,ikbu) ) 
     2775               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp) 
     2776               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp) 
     2777               zwx(ji,jj,ikbu) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,ikbu) 
     2778               ! 
     2779               zfp_vj = vin(ji,jj,ikbv) + ABS( vin(ji,jj,ikbv) ) 
     2780               zfm_vj = vin(ji,jj,ikbv) - ABS( vin(ji,jj,ikbv) ) 
     2781               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp) 
     2782               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp) 
     2783               zwy(ji,jj,ikbv) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,ikbv) 
     2784            END DO 
     2785         END DO 
     2786      ENDIF 
     2787 
     2788      ! total advective trend 
     2789      DO jk = 1, jpkm1 
     2790         DO jj = 2, jpjm1 
     2791            DO ji = fs_2, fs_jpim1   ! vector opt. 
     2792               zbtr = r1_e12t(ji,jj) 
     2793               ! total intermediate advective trends 
     2794               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )  & 
     2795                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )  ) 
     2796               ! 
     2797               ! update and guess with monotonic sheme 
     2798               pta(ji,jj,jk) =   pta(ji,jj,jk) + ztra 
     2799               zwi(ji,jj,jk) = (fse3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     2800            END DO 
     2801         END DO 
     2802      END DO 
     2803      !                             ! Lateral boundary conditions on zwi  (unchanged sign) 
     2804      CALL lbc_lnk( zwi, 'T', 1. )   
     2805 
     2806      IF ( ln_vvl_dbg ) THEN 
     2807         zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 )  
     2808         IF( lk_mpp )   CALL mpp_min( zmin ) 
     2809         IF( zmin < 0._wp) THEN 
     2810            IF(lwp) CALL ctl_stop('vvl_adv: CFL issue here') 
     2811         ENDIF 
     2812      ENDIF 
     2813 
     2814      ! 3. antidiffusive flux : high order minus low order 
     2815      ! -------------------------------------------------- 
     2816      ! antidiffusive flux on i and j 
     2817      DO jk = 1, jpkm1 
     2818         DO jj = 1, jpjm1 
     2819            DO ji = 1, fs_jpim1   ! vector opt. 
     2820               zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * fse3u_n(ji,jj,jk) & 
     2821                                & - zwx(ji,jj,jk)) * umask(ji,jj,jk) 
     2822               zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * fse3v_n(ji,jj,jk) &  
     2823                                & - zwy(ji,jj,jk)) * vmask(ji,jj,jk) 
     2824               ! 
     2825               ! Update advective fluxes 
     2826               un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk) 
     2827               vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk) 
     2828            END DO 
     2829         END DO 
     2830      END DO 
     2831       
     2832      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions 
     2833 
     2834      ! 4. monotonicity algorithm 
     2835      ! ------------------------- 
     2836      CALL nonosc_2d( fse3t_b(:,:,:), zwx, zwy, zwi, z2dtt ) 
     2837 
     2838      ! 5. final trend with corrected fluxes 
     2839      ! ------------------------------------ 
     2840      ! 
     2841      ! Update advective fluxes 
     2842      un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:) 
     2843      vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:) 
     2844      ! 
     2845      DO jk = 1, jpkm1 
     2846         DO jj = 2, jpjm1 
     2847            DO ji = fs_2, fs_jpim1   ! vector opt.   
     2848               ! 
     2849               zbtr = r1_e12t(ji,jj) 
     2850               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     2851                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
     2852               ! add them to the general tracer trends 
     2853               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 
     2854            END DO 
     2855         END DO 
     2856      END DO 
     2857 
     2858      ! 
     2859      ! 6. Correct barotropic flux 
     2860      ! -------------------------- 
     2861      ! Compute barotropic flux difference: 
     2862      zbu(:,:) = 0.e0 
     2863      zbv(:,:) = 0.e0 
     2864      DO jj = 1, jpj 
     2865         DO ji = 1, jpi   ! vector opt. 
     2866            DO jk = 1, jpkm1 
     2867               zbu(ji,jj) = zbu(ji,jj) - un_td(ji,jj,jk) * umask(ji,jj,jk) 
     2868               zbv(ji,jj) = zbv(ji,jj) - vn_td(ji,jj,jk) * vmask(ji,jj,jk) 
     2869            END DO 
     2870         END DO 
     2871      ENDDO 
     2872 
     2873      ! Compute upstream depths: 
     2874      zhu_b(:,:) = 0.e0 
     2875      zhv_b(:,:) = 0.e0 
     2876      IF ( ln_sco ) THEN; ibot=0 ; ELSE ; ibot=1 ; ENDIF 
     2877 
     2878      DO jj = 1, jpjm1 
     2879         DO ji = 1, jpim1   ! vector opt. 
     2880            zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 
     2881            ikbu = mbku(ji,jj) - ibot  
     2882            DO jk = 1, ikbu 
     2883               zfp_hi = fse3t_b(ji  ,jj  ,jk) 
     2884               zfm_hi = fse3t_b(ji+1,jj  ,jk) 
     2885               zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi & 
     2886                            &                  + (1._wp-zdi) * zfm_hi & 
     2887                            &                ) * umask(ji,jj,jk) 
     2888            END DO 
     2889            ! 
     2890            zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 
     2891            ikbv = mbkv(ji,jj) - ibot 
     2892            DO jk = 1, ikbv 
     2893               zfp_hj = fse3t_b(ji  ,jj  ,jk) 
     2894               zfm_hj = fse3t_b(ji  ,jj+1,jk) 
     2895               zhv_b(ji,jj) = zhv_b(ji,jj) + (          zdj  * zfp_hj & 
     2896                            &                  + (1._wp-zdj) * zfm_hj & 
     2897                            &                ) * vmask(ji,jj,jk) 
     2898            END DO 
     2899         END DO 
     2900      END DO 
     2901 
     2902      IF ( .NOT.ln_sco ) THEN 
     2903         ! Correct bottom value 
     2904         ! considering "shelf horizon depth" 
     2905         DO jj = 1, jpjm1 
     2906            DO ji = 1, jpim1   ! vector opt. 
     2907               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 
     2908               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 
     2909               ikbu = mbku(ji,jj) 
     2910               ikbv = mbkv(ji,jj) 
     2911               zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp) 
     2912               zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp) 
     2913               zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp) 
     2914               zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp) 
     2915               zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi & 
     2916                            &                 + (1._wp-zdi) * zfm_hi & 
     2917                            &                ) * umask(ji,jj,ikbu) 
     2918               zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj & 
     2919                            &                 + (1._wp-zdj) * zfm_hj & 
     2920                            &                ) * vmask(ji,jj,ikbv) 
     2921            END DO 
     2922         END DO 
     2923      ENDIF 
     2924 
     2925      ! Corrective barotropic velocity (times hor. scale factor) 
     2926      zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1)) 
     2927      zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1)) 
     2928 
     2929      CALL lbc_lnk( zbu(:,:), 'U', -1. ) 
     2930      CALL lbc_lnk( zbv(:,:), 'V', -1. ) 
     2931      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions 
     2932       
     2933      ! Set corrective fluxes in upstream direction: 
     2934      ! 
     2935      zwx(:,:,:) = 0.e0 
     2936      zwy(:,:,:) = 0.e0 
     2937      DO jj = 1, jpjm1 
     2938         DO ji = 1, fs_jpim1   ! vector opt. 
     2939            ! upstream scheme 
     2940            zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 
     2941            zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 
     2942            zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 
     2943            zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 
     2944            DO jk = 1, jpkm1 
     2945               zfp_hi = fse3t_b(ji  ,jj  ,jk) 
     2946               zfm_hi = fse3t_b(ji+1,jj  ,jk) 
     2947               ! 
     2948               zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2949 
     2950               zfp_hj = fse3t_b(ji  ,jj  ,jk) 
     2951               zfm_hj = fse3t_b(ji  ,jj+1,jk) 
     2952               ! 
     2953               zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2954            END DO 
     2955         END DO 
     2956      END DO 
     2957 
     2958      IF ( .NOT.ln_sco ) THEN 
     2959      ! Bottom correction: 
     2960      DO jj = 1, jpjm1 
     2961         DO ji = 1, fs_jpim1   ! vector opt. 
     2962            ! upstream scheme 
     2963            ikbu = mbku(ji,jj) 
     2964            ikbv = mbkv(ji,jj) 
     2965            ! 
     2966            zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 
     2967            zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 
     2968            zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 
     2969            zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 
     2970            ! 
     2971            zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbu), 0._wp) 
     2972            zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,ikbu), 0._wp) 
     2973            ! 
     2974            zwx(ji,jj,ikbu) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) 
     2975            ! 
     2976            zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,ikbv), 0._wp) 
     2977            zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,ikbv), 0._wp) 
     2978            ! 
     2979            zwy(ji,jj,ikbv) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) 
     2980         END DO 
     2981      END DO 
     2982      ENDIF 
     2983 
     2984      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions 
     2985 
     2986      un_td(:,:,:) = un_td(:,:,:) + zwx(:,:,:) 
     2987      vn_td(:,:,:) = vn_td(:,:,:) + zwy(:,:,:) 
     2988      ! 
     2989      ! Update trend with corrective fluxes: 
     2990      DO jk = 1, jpkm1 
     2991         DO jj = 2, jpjm1 
     2992            DO ji = fs_2, fs_jpim1   ! vector opt.   
     2993               ! 
     2994               zbtr = r1_e12t(ji,jj) 
     2995               ! total advective trends 
     2996               ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     2997                  &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
     2998               ! add them to the general tracer trends 
     2999               pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 
     3000            END DO 
     3001         END DO 
     3002      END DO 
     3003 
     3004      CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv) 
     3005      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi) 
     3006      ! 
     3007      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_adv_fct') 
     3008      ! 
     3009   END SUBROUTINE dom_vvl_adv_fct 
     3010 
     3011   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt ) 
     3012      !!--------------------------------------------------------------------- 
     3013      !!                    ***  ROUTINE nonosc_2d  *** 
     3014      !!      
     3015      !! **  Purpose :   compute monotonic thickness fluxes from the upstream  
     3016      !!       scheme and the before field by a nonoscillatory algorithm  
     3017      !! 
     3018      !! **  Method  :   ... ??? 
     3019      !!       warning : pbef and paft must be masked, but the boundaries 
     3020      !!       conditions on the fluxes are not necessary zalezak (1979) 
     3021      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
     3022      !!       in-space based differencing for fluid 
     3023      !!---------------------------------------------------------------------- 
     3024      ! 
     3025      !!---------------------------------------------------------------------- 
     3026      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
     3027      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field 
     3028      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb        ! monotonic fluxes in the 3 directions 
     3029      ! 
     3030      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     3031      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars 
     3032      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
     3033      REAL(wp) ::   zupip1, zupim1, zupjp1, zupjm1, zupb, zupa 
     3034      REAL(wp) ::   zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa 
     3035      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo 
     3036      !!---------------------------------------------------------------------- 
     3037      ! 
     3038      IF( nn_timing == 1 )  CALL timing_start('nonosc2') 
     3039      ! 
     3040      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 
     3041      ! 
     3042 
     3043      zbig  = 1.e+40_wp 
     3044      zrtrn = 1.e-15_wp 
     3045      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp 
     3046 
     3047 
     3048      ! Search local extrema 
     3049      ! -------------------- 
     3050      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
     3051      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   & 
     3052         &        paft * tmask - zbig * ( 1.e0 - tmask )  ) 
     3053      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   & 
     3054         &        paft * tmask + zbig * ( 1.e0 - tmask )  ) 
     3055 
     3056      DO jk = 1, jpkm1 
     3057         z2dtt = p2dt 
     3058         DO jj = 2, jpjm1 
     3059            DO ji = fs_2, fs_jpim1   ! vector opt. 
     3060 
     3061               ! search maximum in neighbourhood 
     3062               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
     3063                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   & 
     3064                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  )) 
     3065 
     3066               ! search minimum in neighbourhood 
     3067               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   & 
     3068                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   & 
     3069                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  )) 
     3070 
     3071               ! positive part of the flux 
     3072               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   & 
     3073                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )    
     3074 
     3075               ! negative part of the flux 
     3076               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   & 
     3077                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) ) 
     3078 
     3079               ! up & down beta terms 
     3080               zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt 
     3081               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
     3082               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
     3083            END DO 
     3084         END DO 
     3085      END DO 
     3086 
     3087      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
     3088 
     3089      ! 3. monotonic flux in the i & j direction (paa & pbb) 
     3090      ! ---------------------------------------- 
     3091      DO jk = 1, jpkm1 
     3092         DO jj = 2, jpjm1 
     3093            DO ji = fs_2, fs_jpim1   ! vector opt. 
     3094               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
     3095               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
     3096               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) ) 
     3097               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu ) 
     3098 
     3099               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
     3100               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
     3101               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 
     3102               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv ) 
     3103            END DO 
     3104         END DO 
     3105      END DO 
     3106      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign) 
     3107      ! 
     3108      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 
     3109      ! 
     3110      IF( nn_timing == 1 )  CALL timing_stop('nonosc2') 
     3111      ! 
     3112   END SUBROUTINE nonosc_2d 
     3113 
    14033114   !!====================================================================== 
    14043115END MODULE domvvl 
    1405  
    1406  
    1407  
Note: See TracChangeset for help on using the changeset viewer.