Changeset 9529


Ignore:
Timestamp:
2018-04-30T16:38:28+02:00 (2 years ago)
Author:
jchanut
Message:

ztilde stability update

Location:
branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/SHARED/field_def.xml

    r9353 r9529  
    9696         <field id="tpt_dep"      long_name="T-point depth"                  standard_name="depth_below_geoid"   unit="m"   grid_ref="grid_T_3D" /> 
    9797         <field id="e3tdef"       long_name="T-cell thickness deformation"                                       unit="%"   grid_ref="grid_T_3D" /> 
     98         <field id="gdepwt"       long_name="T-cell interface depth"                                             unit="m"   grid_ref="grid_T_3D" /> 
    9899         <field id="stiff_tilde"  long_name="Stiffness"                                                          unit=" "   grid_ref="grid_T_3D" /> 
    99100         <field id="depw_tilde"   long_name="Interface displacement"                                             unit="m"   grid_ref="grid_W_3D" /> 
     
    408409         <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_U_3D" /> 
    409410         <field id="un_lf_tilde"  long_name="Sea Water X Velocity (low pass)"                        standard_name="low_pass_sea_water_x_velocity"     unit="m/s"        grid_ref="grid_U_3D" /> 
    410  
    411411 
    412412         <!-- variables available with MLE --> 
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/CONFIG/cfg.txt

    r9353 r9529  
    1212ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1313ORCA2LIM3_16 OPA_SRC LIM_SRC_3 NST_SRC 
     14ORCA2LIM3_LONG OPA_SRC LIM_SRC_3 NST_SRC 
     15RIDGE OPA_SRC 
    1416COMODO_IW OPA_SRC 
    15 ORCA2LIM3_LONG OPA_SRC LIM_SRC_3 NST_SRC 
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r8473 r9529  
    157157         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    158158      CALL iom_put("tpt_dep", fsdept_n(:,:,:) ) 
    159  
    160  
     159      IF( iom_use("gdepwt") )  THEN 
     160         z3d(:,:,1) = (gdepw_n(:,:,1)-sshn(:,:))*tmask(:,:,1) 
     161         DO jk=2,jpk 
     162            z3d(:,:,jk) = (gdepw_n(:,:,jk)-sshn(:,:))*tmask(:,:,jk-1) 
     163         END DO 
     164         CALL iom_put("gdepwt" , z3d(:,:,:) ) 
     165      END IF 
    161166 
    162167      CALL iom_put( "ssh" , sshn )                 ! sea surface height 
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r9353 r9529  
    6060   LOGICAL                                               :: ln_vvl_blp                          ! Bilaplacian thickness diffusion  
    6161   LOGICAL                                               :: ln_vvl_regrid                       ! ensure layer separation  
     62   LOGICAL                                               :: ll_shorizd=.FALSE.                  ! Use "shelf horizon depths"          
    6263   !                                                                                            ! conservation: not used yet 
    6364   REAL(wp)                                              :: rn_ahe3_lap               ! thickness diffusion coefficient (Laplacian) 
     
    6667   REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days] 
    6768   REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation 
     69   REAL(wp)                                              :: hsmall=0.01               ! small thickness 
    6870 
    6971   !! * Module variables 
    7072   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          
     73   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_lf, vn_lf, hdivn_lf    ! 1st order filtered arrays    
    7274   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n  ! baroclinic scale factors 
    7375   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 
     
    226228         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp ) 
    227229         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
    228          ! tendency mask: 
    229230         ! 
    230231         IF( ln_vvl_ztilde_as_zstar ) THEN 
     
    331332      INTEGER                                :: ib, ib_bdy, ip, jp    !   "     "     " 
    332333      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers 
     334      INTEGER                                :: ncall 
    333335      REAL(wp)                               :: z2dt                  ! temporary scalars 
    334336      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars 
    335337      REAL(wp)                               :: zalpha, zwgt          ! temporary scalars 
    336       REAL(wp)                               :: zdu, zdv 
     338      REAL(wp)                               :: zdu, zdv, zramp 
    337339      LOGICAL                                :: ll_do_bclinic         ! temporary logical 
     340      REAL(wp)                               :: ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv 
    338341      !!---------------------------------------------------------------------- 
    339342      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt') 
     
    348351 
    349352      ll_do_bclinic = .TRUE. 
     353      ncall=1 
    350354      IF( PRESENT(kcall) ) THEN 
    351          IF ( kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE. 
     355         ! comment line below if tilda coordinate has to be computed at each call 
     356         IF (kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE. 
     357         IF (kcall == 2) ncall=2     
    352358      ENDIF 
    353359 
     360      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     361         z2dt =  rdt 
     362      ELSE 
     363         z2dt = 2.0_wp * rdt 
     364      ENDIF 
     365 
    354366      ! ******************************* ! 
    355       ! After acale factors at t-points ! 
     367      ! After scale factors at t-points ! 
    356368      ! ******************************* ! 
    357369 
     
    383395         IF ( ln_vvl_ztilde ) THEN 
    384396            ! 
    385             IF ((kt==nit000).AND.(neuler==0)) THEN 
    386                DO jk = 1, jpkm1 
    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) 
    395                END DO 
    396             ENDIF 
     397            DO jk = 1, jpkm1 
     398               ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk) 
     399               ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk) 
     400               tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
     401            END DO 
     402 
    397403            ! 
    398404         ELSEIF ( ln_vvl_layer ) THEN 
     
    405411         ENDIF 
    406412         ! 
     413         ! Block fluxes through small layers: 
     414!!         DO jk=1,jpkm1 
     415!!            DO ji = 1, jpi 
     416!!               DO jj= 1, jpj 
     417!!                  zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, fse3u_n(ji,jj,jk) - hsmall) ) 
     418!!                  un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * fse3u_n(ji,jj,jk) * e2u(ji,jj) 
     419!!                  ztu(ji,jj,jk) = zmsku * ztu(ji,jj,jk) 
     420!!                  IF ( ln_vvl_ztilde ) un_lf(ji,jj,jk) = zmsku * un_lf(ji,jj,jk) 
     421!!                  ! 
     422!!                  zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, fse3v_n(ji,jj,jk) - hsmall) ) 
     423!!                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * fse3v_n(ji,jj,jk) * e1v(ji,jj) 
     424!!                  ztv(ji,jj,jk) = zmskv * ztv(ji,jj,jk) 
     425!!                  IF ( ln_vvl_ztilde ) vn_lf(ji,jj,jk) = zmskv * vn_lf(ji,jj,jk) 
     426!!               END DO 
     427!!            END DO 
     428!!         END DO       
     429         ! 
    407430         ! Do advection 
    408431         IF     (ln_vvl_adv_fct) THEN 
     
    411434         ELSEIF (ln_vvl_adv_cn2) THEN 
    412435            DO jk = 1, jpkm1 
    413                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * hdivn(:,:,jk) 
     436               DO jj = 2, jpjm1 
     437                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     438                     tilde_e3t_a(ji,jj,jk) =   & 
     439                    -(  e2u(ji,jj)*fse3u(ji,jj,jk) * ztu(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * ztu(ji-1,jj  ,jk)       & 
     440                      + e1v(ji,jj)*fse3v(ji,jj,jk) * ztv(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * ztv(ji  ,jj-1,jk)  )    & 
     441                     / ( e1t(ji,jj) * e2t(ji,jj) ) 
     442                  END DO 
     443               END DO 
    414444            END DO 
    415445         ENDIF 
    416446         !  
    417          ! Thickness anomlaly diffusion: 
     447         ! Thickness anomaly diffusion: 
    418448         ! ----------------------------- 
    419          zwu(:,:) = 0.0_wp 
    420          zwv(:,:) = 0.0_wp 
    421449         ztu(:,:,:) = 0.0_wp 
    422450         ztv(:,:,:) = 0.0_wp 
     451 
     452         ze3t(:,:,1) = 0.e0 
     453         DO jj = 1, jpj 
     454            DO ji = 1, jpi 
     455               DO jk = 2, jpk 
     456                  ze3t(ji,jj,jk) =  ze3t(ji,jj,jk-1) + tilde_e3t_b(ji,jj,jk-1) * tmask(ji,jj,jk-1)  
     457               END DO 
     458            END DO 
     459         END DO 
    423460 
    424461         IF ( ln_vvl_blp ) THEN  ! Bilaplacian 
     
    426463               DO jj = 1, jpjm1                 ! First derivative (gradient) 
    427464                  DO ji = 1, fs_jpim1   ! vector opt.                   
     465!                     ztu(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     466!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     467!                     ztv(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     468!                                     &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    428469                     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) ) 
     470                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) ) 
    430471                     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) ) 
     472                                     &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) ) 
    432473                  END DO 
    433474               END DO 
     
    455496                     ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( zht(ji+1,jj  ) - zht(ji,jj) ) 
    456497                     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) 
    459498                  END DO 
    460499               END DO 
     
    466505               DO jj = 1, jpjm1 
    467506                  DO ji = 1, fs_jpim1   ! vector opt.                   
     507!                     zdu =  rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     508!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     509!                     zdv =  rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     510!                         &  * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    468511                     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) ) 
     512                         &  * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj  ,jk) ) 
    470513                     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 
     514                         &  * ( ze3t(ji,jj,jk) - ze3t(ji  ,jj+1,jk) ) 
    474515                     ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu 
    475516                     ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv 
     
    478519            END DO 
    479520         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 
    496521 
    497522         ! divergence of diffusive fluxes 
     
    499524            DO jj = 2, jpjm1 
    500525               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)    & 
     526!                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk) - ztu(ji,jj,jk)    & 
     527!                     &                                          +     ztv(ji  ,jj-1,jk) - ztv(ji,jj,jk)    & 
     528!                     &                                            ) * r1_e12t(ji,jj) 
     529                  un_td(ji,jj,jk) = un_td(ji,jj,jk) + ztu(ji,jj,jk+1) - ztu(ji,jj,jk  )  
     530                  vn_td(ji,jj,jk) = vn_td(ji,jj,jk) + ztv(ji,jj,jk+1) - ztv(ji,jj,jk  )  
     531                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   ztu(ji-1,jj  ,jk+1) - ztu(ji,jj,jk+1)    & 
     532                     &                                               +ztv(ji  ,jj-1,jk+1) - ztv(ji,jj,jk+1)    & 
     533                     &                                               -ztu(ji-1,jj  ,jk  ) + ztu(ji,jj,jk  )    & 
     534                     &                                               -ztv(ji  ,jj-1,jk  ) + ztv(ji,jj,jk  )    & 
    503535                     &                                            ) * r1_e12t(ji,jj) 
    504536               END DO 
    505537            END DO 
    506538         END DO 
     539 
    507540  
    508          un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:) 
    509          vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:) 
     541!         un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:) 
     542!         vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:) 
     543 
    510544         CALL lbc_lnk( un_td , 'U' , -1.) 
    511          CALL lbc_lnk( vn_td , 'V' , -1.)   
    512          ! 
    513          ! 
    514          ! Restoring: 
    515          IF( ln_vvl_ztilde ) THEN 
    516             DO jk = 1, jpk 
    517                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
    518             END DO 
    519          ENDIF 
    520  
     545         CALL lbc_lnk( vn_td , 'V' , -1.)  
     546         ! 
     547         CALL dom_vvl_ups_cor( kt, tilde_e3t_a, un_td, vn_td ) 
     548 
     549!         IF ( ln_vvl_ztilde ) THEN 
     550!            ztu(:,:,:) = -un_lf(:,:,:) 
     551!            ztv(:,:,:) = -vn_lf(:,:,:) 
     552!            CALL dom_vvl_ups_cor( kt, tilde_e3t_a, ztu, ztv ) 
     553!         ENDIF 
     554         ! 
    521555         ! Remove "external thickness" tendency: 
    522556         DO jk = 1, jpkm1 
     
    526560         ! Leapfrog time stepping 
    527561         ! ~~~~~~~~~~~~~~~~~~~~~~ 
    528          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    529             z2dt =  rdt 
    530          ELSE 
    531             z2dt = 2.0_wp * rdt 
    532          ENDIF 
    533  
    534          tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
    535  
    536          ! Revert to zstar locally: 
    537          ! ~~~~~~~~~~~~~~~~~~~~~~~~ 
     562         zramp = 1._wp 
     563!         IF (.NOT.ln_rstart) zramp = MIN(MAX( ((kt-nit000)*rdt)/(3._wp*rday),0._wp),1._wp) 
     564 
    538565         DO jk=1,jpkm1 
    539             tilde_e3t_a(:,:,jk) =  tilde_e3t_a(:,:,jk) * tildemask(:,:)  
     566            tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) & 
     567                               & * tildemask(:,:) * zramp 
    540568         END DO 
    541569 
     
    595623               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 
    596624            ENDIF 
    597          ENDIF 
    598  
    599          IF ( ln_vvl_ztilde ) THEN 
     625         ENDIF          
     626      ENDIF 
     627 
     628      IF( ln_vvl_ztilde )  THEN 
     629         IF ( ncall==1 ) THEN 
    600630            zalpha  = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
    601631            DO jk = 1, jpkm1 
     
    603633               ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk) 
    604634               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  
    615635            END DO 
    616636            ! 
    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(:,:,:) 
     637               un_lf(:,:,:)  =    un_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:) 
     638               vn_lf(:,:,:)  =    vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:) 
     639            hdivn_lf(:,:,:)  = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:) 
     640!               un_lf(:,:,:)  =    un_lf(:,:,:) * (1._wp - zalpha) + zalpha * un_lf2(:,:,:) 
     641!               vn_lf(:,:,:)  =    vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * vn_lf2(:,:,:) 
     642!            hdivn_lf(:,:,:)  = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * hdivn_lf2(:,:,:) 
    620643         ENDIF 
    621  
     644      ENDIF 
     645 
     646      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
     647      !                                             ! ---baroclinic part--------- ! 
     648         DO jk = 1, jpkm1 
     649            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                      
     650         END DO 
    622651      ENDIF 
    623652 
     
    639668                                  & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:)) 
    640669         END DO 
     670         DO jk = 1, jpkm1 
     671            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * &  
     672                                  & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:)) 
     673         END DO 
    641674      ENDIF 
    642675 
    643       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
    644       !                                             ! ---baroclinic part--------- ! 
    645          DO jk = 1, jpkm1 
    646             fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))                      
    647          END DO 
    648       ENDIF 
    649  
    650       IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging 
     676 
     677      IF( ln_vvl_dbg .AND. ( ncall==2 ) ) THEN   ! - ML - test: control prints for debuging 
    651678         ! 
    652679         zht(:,:) = 0.0_wp 
     
    663690         z_tmin = MINVAL( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0  ) 
    664691         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    665          IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin          
     692         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin 
     693         IF ( z_tmin .LE. 0._wp ) THEN 
     694            IF( lk_mpp ) THEN 
     695               CALL mpp_minloc(fse3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     696            ELSE 
     697               ijk_min = MINLOC( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0   ) 
     698               ijk_min(1) = ijk_min(1) + nimpp - 1 
     699               ijk_min(2) = ijk_min(2) + njmpp - 1 
     700            ENDIF 
     701            IF (lwp) THEN 
     702               WRITE(numout, *) 'Negative scale factor, fse3t_n =', z_tmin 
     703               WRITE(numout, *) 'at i, j, k=', ijk_min             
     704               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
     705            ENDIF 
     706         ENDIF          
    666707         ! 
    667708         z_tmin = MINVAL( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 ) 
    668709         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    669710         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3u_n) =', z_tmin 
     711         IF ( z_tmin .LE. 0._wp ) THEN 
     712            IF( lk_mpp ) THEN 
     713               CALL mpp_minloc(fse3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     714            ELSE 
     715               ijk_min = MINLOC( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0   ) 
     716               ijk_min(1) = ijk_min(1) + nimpp - 1 
     717               ijk_min(2) = ijk_min(2) + njmpp - 1 
     718            ENDIF 
     719            IF (lwp) THEN 
     720               WRITE(numout, *) 'Negative scale factor, fse3u_n =', z_tmin 
     721               WRITE(numout, *) 'at i, j, k=', ijk_min             
     722               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
     723            ENDIF 
     724         ENDIF  
    670725         ! 
    671726         z_tmin = MINVAL( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 ) 
    672727         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    673728         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3v_n) =', z_tmin 
     729         IF ( z_tmin .LE. 0._wp ) THEN 
     730            IF( lk_mpp ) THEN 
     731               CALL mpp_minloc(fse3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     732            ELSE 
     733               ijk_min = MINLOC( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0   ) 
     734               ijk_min(1) = ijk_min(1) + nimpp - 1 
     735               ijk_min(2) = ijk_min(2) + njmpp - 1 
     736            ENDIF 
     737            IF (lwp) THEN 
     738               WRITE(numout, *) 'Negative scale factor, fse3v_n =', z_tmin 
     739               WRITE(numout, *) 'at i, j, k=', ijk_min             
     740               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
     741            ENDIF 
     742         ENDIF  
    674743         ! 
    675744         z_tmin = MINVAL( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 ) 
    676745         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    677746         IF( lwp    ) WRITE(numout, *) kt,' MINVAL(fse3f_n) =', z_tmin 
     747         IF ( z_tmin .LE. 0._wp ) THEN 
     748            IF( lk_mpp ) THEN 
     749               CALL mpp_minloc(fse3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
     750            ELSE 
     751               ijk_min = MINLOC( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0   ) 
     752               ijk_min(1) = ijk_min(1) + nimpp - 1 
     753               ijk_min(2) = ijk_min(2) + njmpp - 1 
     754            ENDIF 
     755            IF (lwp) THEN 
     756               WRITE(numout, *) 'Negative scale factor, fse3f_n =', z_tmin 
     757               WRITE(numout, *) 'at i, j, k=', ijk_min             
     758               CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 
     759            ENDIF 
     760         ENDIF 
    678761         ! 
    679762         zht(:,:) = 0.0_wp 
     
    9561039         zout(:,:,jpk)=0._wp 
    9571040         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(:,:,:) 
     1041            DO jj = 1, jpjm1 
     1042               DO ji = 1, fs_jpim1   ! vector opt.                   
     1043                  ua(ji,jj,jk) =  umask(ji,jj,jk) * re2u_e1u(ji,jj) & 
     1044                                  &  * ( zwdw(ji,jj,jk) - zwdw(ji+1,jj  ,jk) ) 
     1045                  va(ji,jj,jk) =  vmask(ji,jj,jk) * re1v_e2v(ji,jj) &  
     1046                                  &  * ( zwdw(ji,jj,jk) - zwdw(ji  ,jj+1,jk) ) 
     1047               END DO 
     1048            END DO 
     1049         END DO 
     1050             
    9811051         DO jk = 2, jpkm1 
    9821052            DO jj = 2, jpjm1 
    9831053               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 
     1054                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    & 
     1055                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * SQRT(r1_e12t(ji,jj)) 
     1056                  zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)            
     1057               END DO 
     1058            END DO 
     1059         END DO  
    9941060         ! Mask open boundaries: 
    9951061#if defined key_bdy 
     
    10651131      INTEGER ::   jkbot                                               ! bottom level 
    10661132      LOGICAL ::   l_is_orca                                           ! local logical 
    1067       REAL(wp) :: zmin, zdo, zup, ztap 
     1133      REAL(wp) :: zmin, zdo, zup, ztap, zsmall 
    10681134      REAL(wp), POINTER, DIMENSION(:,:)   :: zs                        ! Surface interface depth anomaly 
    10691135      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw                        ! Interface depth anomaly 
     
    10751141 
    10761142  
    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 
     1143!      nmet=0 ! Original method (Surely wrong) 
     1144!      nmet= 1 ! Interface interpolation 
     1145      nmet=1 ! Internal interfaces interpolation only, spread barotropic increment 
     1146 
     1147      ztap = 0._wp ! Minimum fraction of T-point thickness at cell interfaces 
     1148      zsmall = 1e-3_wp 
    10821149 
    10831150      IF ( (nmet==1).OR.(nmet==2) ) THEN 
     
    10911158            zw(:,:,:) =  0._wp     
    10921159            ! 
    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 
     1160!            DO jj = 1, jpj 
     1161!               DO ji = 1, jpi 
     1162!                  jkbot = mbkt(ji,jj) 
     1163!                  DO jk=jkbot,1,-1 
     1164!                     zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk) 
     1165!                  END DO 
     1166!               END DO 
     1167!            END DO  
     1168            ! 
     1169            DO jk=2,jpk 
     1170               zw(:,:,jk) = zw(:,:,jk-1) + pe3_in(:,:,jk-1)*tmask(:,:,jk-1) 
    11001171            END DO  
    1101             ! 
     1172            ! Interface depth anomalies: 
     1173            DO jk=1,jpkm1 
     1174               zw(:,:,jk) = zw(:,:,jk) - zw(:,:,jpk) + ht_0(:,:) 
     1175            END DO 
     1176            zw(:,:,jpk) = ht_0(:,:) 
    11021177            ! 
    11031178            IF (nmet==2) THEN        ! Consider "internal" interfaces only 
    11041179               zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh) 
    11051180               ! 
    1106                zw(:,:,:) = 0._wp 
    11071181               DO jj = 1, jpj 
    11081182                  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)                     & 
     1183                     DO jk=1,jpk 
     1184                        zw(ji,jj,jk) = (zw(ji,jj,jk) + zs(ji,jj))                                          & 
    11121185                                     & * 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) 
     1186                                     & * tmask(ji,jj,jk) 
    11141187                     END DO 
    11151188                  END DO 
     
    11441217                  ! Correction at last level: 
    11451218                  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)) 
     1219                  zdo   = hu_0(ji,jj) 
     1220                  DO jk=jkbot,1,-1 
     1221                     zup = 0.5_wp * umask(ji,jj,jk)   * r1_e12u(ji,jj)   & 
     1222                           &      * (   e12t(ji  ,jj) * zw(ji  ,jj,jk)   & 
     1223                           &          + e12t(ji+1,jj) * zw(ji+1,jj,jk)   ) 
     1224                     ! 
     1225                     ! If there is a step, taper bottom interface: 
     1226                     IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(zup>zdo)) THEN 
     1227                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN 
     1228!                           zmin = ztap * pe3_in(ji+1,jj,jk) 
     1229                           zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk)) 
     1230                        ELSE 
     1231!                           zmin = ztap * pe3_in(ji  ,jj,jk) 
     1232                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk)) 
     1233                        ENDIF 
     1234                        zup = MIN(zup, zdo-zmin) 
    11581235                     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) 
     1236                     zup = MIN(zup, zdo-zsmall) 
     1237                     pe3_out(ji,jj,jk) = zdo - zup - e3u_0(ji,jj,jk) 
     1238                     zdo = zup 
    11701239                  END DO 
    11711240               END DO 
     
    12171286                  ! Correction at last level: 
    12181287                  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)) 
     1288                  zdo   = hv_0(ji,jj) 
     1289                  DO jk=jkbot,1,-1 
     1290                     zup = 0.5_wp * vmask(ji,jj,jk)   * r1_e12v(ji,jj)  & 
     1291                           &      * (   e12t(ji,jj  ) * zw(ji,jj  ,jk)   & 
     1292                           &          + e12t(ji,jj+1) * zw(ji,jj+1,jk)   ) 
     1293                     ! 
     1294                     ! If there is a step, taper bottom interface: 
     1295                     IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN 
     1296                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN 
     1297!                           zmin = ztap * pe3_in(ji,jj+1,jk) 
     1298                           zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk)) 
     1299                        ELSE 
     1300!                           zmin = ztap * pe3_in(ji  ,jj,jk) 
     1301                           zmin = ztap * (zw(ji,jj  ,jk+1)-zw(ji,jj  ,jk)) 
     1302                        ENDIF 
     1303                        zup = MIN(zup, zdo-zmin)                         
    12311304                     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) 
     1305                     zup = MIN(zup, zdo-zsmall) 
     1306                     pe3_out(ji,jj,jk) = zdo - zup - e3v_0(ji,jj,jk) 
     1307                     zdo = zup 
    12441308                  END DO 
    12451309               END DO 
     
    12941358                  ! bottom correction: 
    12951359                  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 
     1360                  zdo = hf_0(ji,jj) 
     1361                  DO jk=jkbot,1,-1 
    13271362                     zup =  0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 
    13281363                           &        * (  e12t(ji  ,jj  ) * zw(ji  ,jj  ,jk)  &  
     
    13301365                           &           + e12t(ji  ,jj+1) * zw(ji  ,jj+1,jk)  &  
    13311366                           &           + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jk)  ) 
    1332                      pe3_out(ji,jj,jk) = zdo - zup  
    13331367                     ! 
    1334                      zdo =  zdo - pe3_out(ji,jj,jk)                                     
     1368                     ! If there is a step, taper bottom interface: 
     1369                     IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN 
     1370                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN 
     1371                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN 
     1372!                              zmin = ztap * pe3_in(ji+1,jj+1,jk) 
     1373                              zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk)) 
     1374                           ELSE 
     1375!                              zmin = ztap * pe3_in(ji  ,jj+1,jk) 
     1376                              zmin = ztap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk)) 
     1377                           ENDIF 
     1378                        ELSE 
     1379                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN 
     1380!                              zmin = ztap * pe3_in(ji+1,jj  ,jk) 
     1381                              zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk)) 
     1382                           ELSE 
     1383!                              zmin = ztap * pe3_in(ji  ,jj  ,jk) 
     1384                              zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk)) 
     1385                           ENDIF 
     1386                        ENDIF 
     1387                        zup = MIN(zup, zdo-zmin) 
     1388                     ENDIF 
     1389                     zup = MIN(zup, zdo-zsmall) 
     1390                     ! 
     1391                     pe3_out(ji,jj,jk) = zdo - zup - e3f_0(ji,jj,jk) 
     1392                     zdo = zup 
    13351393                  END DO 
    1336                   ! 
    1337                END DO 
    1338             END DO  
     1394               END DO 
     1395            END DO 
    13391396         ENDIF 
    13401397         ! 
     
    14321489            id(3) = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
    14331490            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  
     1491            id(5) = iom_varid( numror, 'un_lf'   , ldstop = .FALSE. ) 
     1492            id(6) = iom_varid( numror, 'vn_lf'   , ldstop = .FALSE. ) 
     1493            id(7) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. ) 
    14381494            !                             ! --------- ! 
    14391495            !                             ! all cases ! 
     
    14981554                  !                       ! ------------ ! 
    14991555                  IF( MINVAL(id(5:7) ) > 0 ) THEN  ! all required arrays exist 
    1500                      CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) ) 
    15011556                     CALL iom_get( numror, jpdom_autoglo, 'un_lf'   ,    un_lf(:,:,:) ) 
    15021557                     CALL iom_get( numror, jpdom_autoglo, 'vn_lf'   ,    vn_lf(:,:,:) ) 
     1558                     CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) ) 
    15031559                  ELSE                            ! one at least array is missing 
    1504                      hdivn_lf(:,:,:) = 0.0_wp 
    15051560                     un_lf(:,:,:)    = 0.0_wp 
    15061561                     vn_lf(:,:,:)    = 0.0_wp 
     1562                     hdivn_lf(:,:,:) = 0.0_wp 
    15071563                     neuler = 0 
    15081564                  ENDIF 
     
    15201576            END IF 
    15211577            IF( ln_vvl_ztilde ) THEN 
    1522                hdivn_lf(:,:,:) = 0.0_wp 
    15231578               un_lf(:,:,:)    = 0.0_wp 
    15241579               vn_lf(:,:,:)    = 0.0_wp 
     1580               hdivn_lf(:,:,:) = 0.0_wp 
    15251581            ENDIF 
    15261582         ENDIF 
     
    15431599         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    15441600            !                                        ! ------------ ! 
    1545             CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:) ) 
    15461601            CALL iom_rstput( kt, nitrst, numrow, 'un_lf'   , un_lf(:,:,:)    ) 
    15471602            CALL iom_rstput( kt, nitrst, numrow, 'vn_lf'   , vn_lf(:,:,:)    ) 
     
    16511706      ENDIF 
    16521707 
     1708      ! Use of "shelf horizon depths" should be allowed with s-z coordinates, but we restrict it to zco and zps  
     1709      ! for the time being 
     1710      IF ( ln_sco ) THEN  
     1711        ll_shorizd=.FALSE. 
     1712      ELSE 
     1713        ll_shorizd=.TRUE. 
     1714      ENDIF 
     1715 
    16531716#if defined key_agrif 
    16541717      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' ) 
     
    20992162                za3 = 0.5_wp * (za1 + za2)  
    21002163                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 
     2164                IF (zdiff>=0.8) THEN 
     2165                   zwt(ji,jj,jk) =  zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk) 
     2166!                   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) 
     2167 
     2168                ELSE 
     2169                   zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk) 
     2170                ENDIF 
    21082171            END DO 
    21092172         END DO 
     
    23002363      !! 
    23012364      !! ** Method  :  More or less adapted from references below. 
    2302       !! 
     2365      !!regrid 
    23032366      !! ** Action  :  Ensure that thickness are above a given value, spaced enough 
    23042367      !!               and revert to Eulerian coordinates near the bottom.       
     
    23182381      INTEGER                             :: ji, jj, jk ! dummy loop indices 
    23192382      LOGICAL                             :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond 
    2320       LOGICAL                             :: ll_e3tsurf_const, ll_zdiff_cond, ll_blpdiff_cond 
     2383      LOGICAL                             :: ll_zdiff_cond, ll_blpdiff_cond 
    23212384      INTEGER                             :: jkbot 
    2322       REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd, dzmin_vvl 
    2323       REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj 
     2385      REAL(wp)                            :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd 
     2386      REAL(wp)                            :: zufim1, zufi, zvfjm1, zvfj, dzmin_int, dzmin_surf 
    23242387      REAL(wp)                            :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn 
    23252388      REAL(wp)                            :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot 
     
    23362399      ll_chk_bot2top   = .TRUE. 
    23372400      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 
     2401      dzmin_int  = 0.1_wp  ! Absolute minimum depth in the interior (in meters) 
     2402      dzmin_surf = 1.0_wp   ! Absolute minimum depth at the surface (in meters) 
     2403      zfrch_stp  = 0.1_wp   ! Maximum fractionnal thickness change in one time step (<= 1.) 
     2404      zfrch_rel  = 0.4_wp   ! Maximum relative thickness change in the vertical (<= 1.) 
     2405      zfrac_bot  = 0.05_wp  ! Fraction of bottom level allowed to change 
     2406      zscal_bot  = 2.0_wp   ! Depth lengthscale 
     2407      ll_zdiff_cond    = .TRUE.  ! Conditionnal vertical diffusion of interfaces 
     2408         zvdiff  = 0.1_wp   ! m 
     2409         zvlim   = 0.5_wp   ! max d2h/dh 
     2410      ll_lapdiff_cond  = .TRUE.  ! Conditionnal Laplacian diffusion of interfaces 
     2411         zhdiff  = 0.2_wp   ! ad. 
     2412         zhlim   = 0.03_wp  ! ad. max lap(z)*e1 
     2413      ll_blpdiff_cond  = .TRUE.  ! Conditionnal Bilaplacian diffusion of interfaces 
     2414         zhdiff2 = 0.2_wp   ! ad. 
     2415!         zhlim2  = 3.e-11_wp  !  max bilap(z) 
     2416         zhlim2  = 0.03_wp  ! ad. max bilap(z)*e1**3 
    23532417      ! ---------------------------------------------------------------------------------------  
    23542418      ! 
     
    23722436                          &     ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj  ,1), & 
    23732437                          &     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 )  
     2438                zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall ) 
    23752439            END DO 
    23762440         END DO 
     
    23852449                  &                           + 4._wp*  zdw(ji  ,jj  )                       ) 
    23862450            END DO 
    2387          END DO 
     2451         END DO          
    23882452 
    23892453         CALL lbc_lnk( dsm(:,:), 'T', 1. ) 
     
    23952459                  jk = i_int_bot(ji,jj) 
    23962460                  hsm(ji,jj) = zfrac_bot * e3w_1d(jk) 
     2461!                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj)) 
    23972462               END DO 
    23982463            END DO 
     
    24022467                  jk = i_int_bot(ji,jj) 
    24032468                  hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk) 
     2469!                  dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj)) 
    24042470               END DO 
    24052471            END DO 
     
    25152581                  ztmp1 = ( (ua(ji-1,jj  ,jk) - ua(ji,jj,jk))    & 
    25162582                     &  +   (va(ji  ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj) 
    2517                   zh2 = MAX(abs(ztmp1)-zhlim2, 0._wp) 
     2583                  zh2 = MAX(abs(ztmp1)-zhlim2*SQRT(r1_e12t(ji,jj))*r1_e12t(ji,jj), 0._wp) 
    25182584                  ztmp = SIGN(zh2, ztmp1) 
    25192585                  zeu2 = zhdiff2 * e12t(ji,jj)*e12t(ji,jj) / 16._wp 
     
    25552621                  ! 
    25562622                  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) 
     2623                  zh_bef = tilde_e3t_b(ji,jj,jk) + zh_0 
     2624                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
     2625                  zh_dwn = tilde_e3t_a(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) 
     2626                  zh_min = MIN(zh_0/3._wp, dzmin_int) 
    25612627                  !  
    25622628                  ! Set maximum and minimum vertical excursions    
    25632629                  ztmph = hsm(ji,jj) 
    25642630                  ztmpd = dsm(ji,jj) 
    2565                   zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)                
     2631                  zh2   = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd) 
    25662632                  zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 ) 
    25672633                  zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff) 
     
    25732639                  ! 
    25742640                  ! Ensure minimum layer thickness:                   
    2575 !                  zh_new = MIN(zh_old, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) ) 
     2641!                  zh_new = MAX(zh_new, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) ) 
     2642                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new) 
    25762643                  zh_new = cush(zh_new, zh_min) 
    25772644                  ! 
     
    25792646                  zdiff = (zh_new - zh_old) * tmask(ji,jj,jk) 
    25802647                  ! 
    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 
     2648                  ! Limit thickness change in 1 time step:  
     2649!                  zh_new = MIN( ABS(zh_new-zh_bef), (1._wp-zfrch_stp)*zh_bef ) 
     2650!                  zdiff = SIGN(ztmp, zh_new - zh_old) 
     2651!                  zh_new = zdiff + zh_old 
    25852652                  ! 
    2586                   tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
     2653!                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
    25872654                  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)             
     2655!                  tilde_e3t_a(ji,jj,jk-1) = (-zdiff + tilde_e3t_a(ji,jj,jk-1) ) * tmask(ji,jj,jk-1)          
    25892656               END DO     
    25902657            END DO  
     
    26012668                  ! 
    26022669                  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) 
     2670                  zh_bef = tilde_e3t_b(ji,jj,jk) + zh_0 
     2671                  zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
     2672                  zh_up  = tilde_e3t_a(ji,jj,jk+1) + e3t_0(ji,jj,jk+1) 
     2673                  zh_min = MIN(zh_0/3._wp, dzmin_int) 
     2674                  ! 
     2675                  zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf) 
     2676                  ! 
     2677                  ! New layer thickness: 
     2678                  zh_new =  zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 
    26082679                  ! 
    26092680                  ! 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) 
     2681!                  zh_new=MAX(zh_new, zh_up * zfrch_rel / (2._wp-zfrch_rel) ) 
     2682                  zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new) 
     2683                  zh_new = cush(zh_new, zh_min) 
    26122684                  ! 
    26132685                  ! Final flux: 
     
    26152687                  !  
    26162688                  ! Limit flux:                  
    2617                   ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 
    2618                   zdiff = SIGN(ztmp, zdiff) 
    2619                   zh_new = zdiff + zh_old 
     2689!                  ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 
     2690!                  zdiff = SIGN(ztmp, zdiff) 
     2691!                  zh_new = zdiff + zh_old 
    26202692                  ! 
    2621                   tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
     2693!                  tilde_e3t_a(ji,jj,jk  ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
    26222694                  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) 
     2695!                  tilde_e3t_a(ji,jj,jk+1) = (-zdiff + tilde_e3t_a(ji,jj,jk+1) ) * tmask(ji,jj,jk+1) 
    26242696               END DO 
    26252697               !                
     
    26272699         END DO 
    26282700      ENDIF 
     2701      ! 
     2702      DO jj = 2, jpjm1 
     2703         DO ji = 2, jpim1 
     2704            DO jk = 1, jpkm1 
     2705               tilde_e3t_a(ji,jj,jk) =  (zwdw(ji,jj,jk+1)-zwdw(ji,jj,jk)-e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 
     2706            END DO 
     2707         END DO 
     2708      END DO 
    26292709      ! 
    26302710      CALL wrk_dealloc( jpi, jpj, jpk, zwdw )  
     
    27032783      !! **  Action  : - Update pta thickness tendency and diffusive fluxes 
    27042784      !!               - 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 
    27072785      !!---------------------------------------------------------------------- 
    27082786      ! 
    27092787      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    27102788      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   
     2789      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input velocities 
     2790      ! 
     2791      INTEGER  ::   ji, jj, jk, ib, ib_bdy               ! dummy loop indices   
    27142792      INTEGER  ::   ikbu, ikbv, ibot 
    27152793      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar 
     
    27192797      REAL(wp) ::   zfp_hi, zfp_hj           !   -      - 
    27202798      REAL(wp) ::   zfm_hi, zfm_hj           !   -      - 
    2721       REAL(wp), POINTER, DIMENSION(:,:)   :: zbu, zbv, zhu_b, zhv_b 
     2799      REAL(wp) ::   ztout,  ztin, zfac       !   -      - 
    27222800      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwi 
    27232801      !!---------------------------------------------------------------------- 
     
    27252803      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_adv_fct') 
    27262804      ! 
    2727       CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv) 
    27282805      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwi) 
    27292806      ! 
     
    27442821      ! 2. upstream advection with initial mass fluxes & intermediate update 
    27452822      ! -------------------------------------------------------------------- 
    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) 
     2823      IF ( ll_shorizd ) THEN 
     2824         DO jk = 1, jpkm1 
     2825            DO jj = 1, jpjm1 
     2826               DO ji = 1, fs_jpim1   ! vector opt. 
     2827                  ! 
     2828                  zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,jk), 0._wp) 
     2829                  zfp_hi = MIN(zfp_hi, fse3t_b(ji  ,jj  ,jk)) 
     2830                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )  
     2831                  ! 
     2832                  zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,jk), 0._wp) 
     2833                  zfm_hi = MIN(zfm_hi, fse3t_b(ji+1,jj  ,jk)) 
     2834                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )  
     2835                  !  
     2836                  zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,jk), 0._wp) 
     2837                  zfp_hj = MIN(zfp_hj, fse3t_b(ji  ,jj  ,jk)) 
     2838                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )  
     2839                  ! 
     2840                  zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,jk), 0._wp) 
     2841                  zfm_hj = MIN(zfm_hj, fse3t_b(ji  ,jj+1,jk)) 
     2842                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 
     2843                  ! 
     2844                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 
     2845                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 
     2846                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 
     2847                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 
     2848                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2849                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2850               END DO 
     2851            END DO 
     2852         END DO 
     2853      ELSE 
     2854         DO jk = 1, jpkm1 
     2855            DO jj = 1, jpjm1 
     2856               DO ji = 1, fs_jpim1   ! vector opt. 
     2857                  ! 
     2858                  zfp_hi = fse3t_b(ji  ,jj  ,jk) 
     2859                  zfm_hi = fse3t_b(ji+1,jj  ,jk)              
     2860                  zfp_hj = fse3t_b(ji  ,jj  ,jk) 
     2861                  zfm_hj = fse3t_b(ji  ,jj+1,jk) 
     2862                  ! 
     2863                  zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 
     2864                  zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 
     2865                  zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 
     2866                  zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 
     2867                  zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     2868                  zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     2869               END DO 
    27842870            END DO 
    27852871         END DO 
     
    28012887         END DO 
    28022888      END DO 
    2803       !                             ! Lateral boundary conditions on zwi  (unchanged sign) 
     2889 
    28042890      CALL lbc_lnk( zwi, 'T', 1. )   
     2891 
     2892#if defined key_bdy 
     2893         DO ib_bdy=1, nb_bdy 
     2894            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 
     2895               ji = idx_bdy(ib_bdy)%nbi(ib,1) 
     2896               jj = idx_bdy(ib_bdy)%nbj(ib,1) 
     2897               DO jk = 1, jpkm1   
     2898                  zwi(ji,jj,jk) = fse3t_a(ji,jj,jk) 
     2899               END DO 
     2900            END DO  
     2901         END DO 
     2902#endif 
    28052903 
    28062904      IF ( ln_vvl_dbg ) THEN 
     
    28082906         IF( lk_mpp )   CALL mpp_min( zmin ) 
    28092907         IF( zmin < 0._wp) THEN 
    2810             IF(lwp) CALL ctl_stop('vvl_adv: CFL issue here') 
     2908            IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here') 
     2909            IF(lwp) WRITE(numout,*) zmin 
    28112910         ENDIF 
    28122911      ENDIF 
     
    28552954         END DO 
    28562955      END DO 
    2857  
    2858       ! 
    2859       ! 6. Correct barotropic flux 
    2860       ! -------------------------- 
     2956      ! 
     2957      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi) 
     2958      ! 
     2959      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_adv_fct') 
     2960      ! 
     2961   END SUBROUTINE dom_vvl_adv_fct 
     2962 
     2963   SUBROUTINE dom_vvl_ups_cor( kt, pta, uin, vin )  
     2964      !!---------------------------------------------------------------------- 
     2965      !!                  ***  ROUTINE dom_vvl_adv_fct  *** 
     2966      !!  
     2967      !! **  Purpose :  Correct for addionnal barotropic fluxes  
     2968      !!                in the upstream direction 
     2969      !! 
     2970      !! **  Method  :    
     2971      !! 
     2972      !! **  Action  : - Update diffusive fluxes uin, vin 
     2973      !!               - Remove divergence from thickness tendency 
     2974      !!---------------------------------------------------------------------- 
     2975      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     2976      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta       ! thickness baroclinic trend  
     2977      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   uin, vin  ! input fluxes 
     2978      INTEGER  ::   ji, jj, jk               ! dummy loop indices   
     2979      INTEGER  ::   ikbu, ikbv, ibot 
     2980      REAL(wp) ::   zbtr, ztra               ! local scalar 
     2981      REAL(wp) ::   zdi, zdj                 !   -      - 
     2982      REAL(wp) ::   zfp_hi, zfp_hj           !   -      - 
     2983      REAL(wp) ::   zfm_hi, zfm_hj           !   -      - 
     2984      REAL(wp) ::   zfp_ui, zfp_vj           !   -      - 
     2985      REAL(wp) ::   zfm_ui, zfm_vj           !   -      - 
     2986      REAL(wp), POINTER, DIMENSION(:,:)   :: zbu, zbv, zhu_b, zhv_b 
     2987      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy 
     2988      !!---------------------------------------------------------------------- 
     2989      ! 
     2990      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_ups_cor') 
     2991      ! 
     2992      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv) 
     2993      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy) 
     2994      ! 
    28612995      ! Compute barotropic flux difference: 
    28622996      zbu(:,:) = 0.e0 
     
    28652999         DO ji = 1, jpi   ! vector opt. 
    28663000            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) 
     3001               zbu(ji,jj) = zbu(ji,jj) - uin(ji,jj,jk) * umask(ji,jj,jk) 
     3002               zbv(ji,jj) = zbv(ji,jj) - vin(ji,jj,jk) * vmask(ji,jj,jk) 
    28693003            END DO 
    28703004         END DO 
     
    28743008      zhu_b(:,:) = 0.e0 
    28753009      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 
     3010 
     3011      IF ( ll_shorizd ) THEN 
    29033012         ! Correct bottom value 
    29043013         ! considering "shelf horizon depth" 
     
    29073016               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 
    29083017               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) 
     3018               DO jk=1, jpkm1 
     3019                  zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,jk), 0._wp) 
     3020                  zfp_hi = MIN(zfp_hi, fse3t_b(ji  ,jj  ,jk)) 
     3021                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )  
     3022                  ! 
     3023                  zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,jk), 0._wp) 
     3024                  zfm_hi = MIN(zfm_hi, fse3t_b(ji+1,jj  ,jk)) 
     3025                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 
     3026                  !  
     3027                  zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,jk), 0._wp) 
     3028                  zfp_hj = MIN(zfp_hj, fse3t_b(ji  ,jj  ,jk)) 
     3029                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )  
     3030                  ! 
     3031                  zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,jk), 0._wp) 
     3032                  zfm_hj = MIN(zfm_hj, fse3t_b(ji  ,jj+1,jk)) 
     3033                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 
     3034                  ! 
     3035                  zhu_b(ji,jj) = zhu_b(ji,jj) + (         zdi  * zfp_hi & 
     3036                               &                 + (1._wp-zdi) * zfm_hi & 
     3037                               &                ) * umask(ji,jj,jk) 
     3038                  zhv_b(ji,jj) = zhv_b(ji,jj) + (         zdj  * zfp_hj & 
     3039                               &                 + (1._wp-zdj) * zfm_hj & 
     3040                               &                ) * vmask(ji,jj,jk) 
     3041               END DO 
     3042            END DO 
     3043         END DO 
     3044      ELSE 
     3045         DO jj = 1, jpjm1 
     3046            DO ji = 1, jpim1   ! vector opt. 
     3047               zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))  
     3048               zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 
     3049               DO jk = 1, jpkm1 
     3050                  zfp_hi = fse3t_b(ji  ,jj  ,jk) 
     3051                  zfm_hi = fse3t_b(ji+1,jj  ,jk) 
     3052                  zfp_hj = fse3t_b(ji  ,jj  ,jk) 
     3053                  zfm_hj = fse3t_b(ji  ,jj+1,jk) 
     3054                  ! 
     3055                  zhu_b(ji,jj) = zhu_b(ji,jj) + (          zdi  * zfp_hi & 
     3056                               &                  + (1._wp-zdi) * zfm_hi & 
     3057                               &                ) * umask(ji,jj,jk) 
     3058                  ! 
     3059                  zhv_b(ji,jj) = zhv_b(ji,jj) + (          zdj  * zfp_hj & 
     3060                               &                  + (1._wp-zdj) * zfm_hj & 
     3061                               &                ) * vmask(ji,jj,jk) 
     3062               END DO 
    29213063            END DO 
    29223064         END DO 
     
    29293071      CALL lbc_lnk( zbu(:,:), 'U', -1. ) 
    29303072      CALL lbc_lnk( zbv(:,:), 'V', -1. ) 
    2931       CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions 
    29323073       
    29333074      ! Set corrective fluxes in upstream direction: 
     
    29353076      zwx(:,:,:) = 0.e0 
    29363077      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 
     3078      IF ( ll_shorizd ) THEN 
     3079         DO jj = 1, jpjm1 
     3080            DO ji = 1, fs_jpim1   ! vector opt. 
     3081               ! upstream scheme 
     3082               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 
     3083               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 
     3084               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 
     3085               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 
     3086               DO jk = 1, jpkm1 
     3087                  zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji  ,jj  ,jk), 0._wp) 
     3088                  zfp_hi = MIN(fse3t_b(ji  ,jj  ,jk), zfp_hi) 
     3089                  zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )  
     3090                  ! 
     3091                  zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj  ,jk), 0._wp) 
     3092                  zfm_hi = MIN(fse3t_b(ji+1,jj  ,jk), zfm_hi) 
     3093                  zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )  
     3094                  ! 
     3095                  zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj  ,jk), 0._wp) 
     3096                  zfp_hj = MIN(fse3t_b(ji  ,jj  ,jk), zfp_hj)  
     3097                  zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )  
     3098                  ! 
     3099                  zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji  ,jj+1,jk), 0._wp) 
     3100                  zfm_hj = MIN(fse3t_b(ji  ,jj+1,jk), zfm_hj) 
     3101                  zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )  
     3102                  ! 
     3103                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     3104                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     3105               END DO 
     3106            END DO 
     3107         END DO 
     3108      ELSE 
     3109         DO jj = 1, jpjm1 
     3110            DO ji = 1, fs_jpim1   ! vector opt. 
     3111               ! upstream scheme 
     3112               zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 
     3113               zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 
     3114               zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 
     3115               zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 
     3116               DO jk = 1, jpkm1 
     3117                  zfp_hi = fse3t_b(ji  ,jj  ,jk) 
     3118                  zfm_hi = fse3t_b(ji+1,jj  ,jk) 
     3119                  zfp_hj = fse3t_b(ji  ,jj  ,jk) 
     3120                  zfm_hj = fse3t_b(ji  ,jj+1,jk) 
     3121                  ! 
     3122                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 
     3123                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 
     3124               END DO 
     3125            END DO 
     3126         END DO 
    29823127      ENDIF 
    2983  
    29843128      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral boundary conditions 
    29853129 
    2986       un_td(:,:,:) = un_td(:,:,:) + zwx(:,:,:) 
    2987       vn_td(:,:,:) = vn_td(:,:,:) + zwy(:,:,:) 
     3130      uin(:,:,:) = uin(:,:,:) + zwx(:,:,:) 
     3131      vin(:,:,:) = vin(:,:,:) + zwy(:,:,:) 
    29883132      ! 
    29893133      ! Update trend with corrective fluxes: 
     
    30013145         END DO 
    30023146      END DO 
    3003  
     3147      ! 
    30043148      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 
     3149      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy) 
     3150      ! 
     3151      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_ups_cor') 
     3152      ! 
     3153   END SUBROUTINE dom_vvl_ups_cor 
    30103154 
    30113155   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt ) 
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r8594 r9529  
    198198      !!---------------------------------------------------------------------- 
    199199      INTEGER  ::   jk                     ! dummy loop indices 
     200      LOGICAL  ::   ll_e3_dep 
    200201      REAL(wp) ::   zt, zw                 ! temporary scalars 
    201202      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in 
     
    209210      ! Set variables from parameters 
    210211      ! ------------------------------ 
     212       ll_e3_dep = .TRUE. 
    211213       zkth = ppkth       ;   zacr = ppacr 
    212214       zdzmin = ppdzmin   ;   zhmax = pphmax 
     
    312314      ENDIF 
    313315 
    314       IF ( ln_isfcav ) THEN 
     316      IF ( ln_isfcav .OR. ll_e3_dep ) THEN 
    315317! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    316318! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r6204 r9529  
    144144 
    145145    
    146    SUBROUTINE wzv( kt ) 
     146   SUBROUTINE wzv( kt, kcall ) 
    147147      !!---------------------------------------------------------------------- 
    148148      !!                ***  ROUTINE wzv  *** 
     
    161161      ! 
    162162      INTEGER, INTENT(in) ::   kt           ! time step 
     163      INTEGER, INTENT( in ), OPTIONAL     :: kcall   ! optional argument 
    163164      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
    164165      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv 
    165166      ! 
     167      LOGICAL             ::   ll_use_totflx 
    166168      INTEGER             ::   ji, jj, jk   ! dummy loop indices 
    167169      REAL(wp)            ::   z1_2dt       ! local scalars 
     
    170172      IF( nn_timing == 1 )  CALL timing_start('wzv') 
    171173      ! 
     174      ll_use_totflx=.FALSE. 
     175      IF (( PRESENT(kcall) ).AND.(ln_vvl_ztilde .OR. ln_vvl_layer)) THEN 
     176         IF ( kcall==2 ) ll_use_totflx=.TRUE. 
     177      ENDIF 
     178 
    172179      IF( kt == nit000 ) THEN 
    173180         ! 
     
    185192      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt 
    186193      ! 
    187       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
     194      IF (( ln_vvl_ztilde .OR. ln_vvl_layer ).AND.ll_use_totflx) THEN      ! z_tilde and layer cases 
    188195         CALL wrk_alloc( jpi, jpj, jpk, zhdiv )  
    189196         ! 
  • branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8477 r9529  
    222222                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 
    223223          IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
    224                                   CALL wzv           ( kstp )  ! now cross-level velocity  
     224                                  CALL wzv           ( kstp, kcall=2 )  ! now cross-level velocity  
    225225      ENDIF 
    226226 
Note: See TracChangeset for help on using the changeset viewer.