Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12492)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12579)
@@ -17,4 +17,5 @@
    !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid
    !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
+   !!   dom_vvl_r3c      : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points
    !!   dom_vvl_rst      : read/write restart file
    !!   dom_vvl_ctl      : Check the vvl options
@@ -245,5 +246,4 @@
       LOGICAL                ::   ll_do_bclinic         ! local logical
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
       !!----------------------------------------------------------------------
       !
@@ -258,8 +258,11 @@
       ENDIF
 
-      ! ll_do_bclinic = .TRUE.
-      ! IF( PRESENT(kcall) ) THEN
-      !    IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
-      ! ENDIF
+      IF( PRESENT(kcall) ) THEN
+         IF( kcall == 2 ) THEN
+            CALL dom_qe_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) )
+         ELSE
+            CALL dom_qe_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa) )
+         ENDIF
+      ENDIF
 
       ! ******************************* !
@@ -274,6 +277,4 @@
       ! After scale factors at u- v- points !
       ! *********************************** !
-      !
-      CALL dom_qe_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa) )
       !
       DO jk = 1, jpkm1
@@ -342,65 +343,46 @@
       ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
       ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
-      !!st  ! r3t/u/v should be unchanged
-      CALL dom_qe_r3c( ssh(:,:,Kmm), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:), r3f(:,:) )
-      !
-      DO jk = 1, jpkm1                    ! Horizontal interpolation of e3t
-         e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:)     * fmask(:,:,jk) )   ! Kmm time level
+
+
+      ! Scale factor computation
+      DO jk = 1, jpk             ! Horizontal interpolation
+         e3f(:,:,jk)      =  e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) )   ! Kmm time level
+         !                       ! Vertical interpolation
+	      !                                   ! The ratio does not have to be masked at w-level
+         e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
+         e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
+         e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
+         e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
+         e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
+         e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
       END DO
-      !CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
-
-      ! Vertical scale factor interpolations
-      ! DO jk = 1, jpk                      ! Vertical interpolation of e3t,u,v
-      !    !                                   ! The ratio does not have to be masked at w-level
-      !    e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
-      !    e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
-      !    e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
-      ! END DO
-      CALL dom_qe_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
-      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
-      CALL dom_qe_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
-      CALL dom_qe_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
-      CALL dom_qe_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
-      CALL dom_qe_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
-
-      ! t- and w- points depth (set the isf depth as it is in the initial step)
-      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
-      gdepw(:,:,1,Kmm) = 0.0_wp
-      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
-      DO_3D_11_11( 2, jpk )
-        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-                                                           ! 1 for jk = mikt
-         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
-         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
-             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )
-         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
-      END_3D
-!       IF( ln_isf ) THEN          !** IceShelF cavities
-!       !                             ! to be created depending of the new names in isf
-!       !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
-!       !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
-! !!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
-!          gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
-!          gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
-!          gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
-!          DO jk = 2, jpk
-!             gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
-!                               + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
-!             gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
-!                               + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
-!             gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
-!          END DO
-!          !
-!       ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
-!          !
-! !!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
-!          DO jk = 1, jpk
-!             gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
-!             gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
-!             gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
-!          END DO
-!          !
-!       ENDIF
+
+
+      IF( ln_isf ) THEN          !** IceShelF cavities
+      !                             ! to be created depending of the new names in isf
+      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
+      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
+!!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
+         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
+         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
+         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+         DO jk = 2, jpk
+            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
+         END DO
+         !
+      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
+         !
+!!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
+         DO jk = 1, jpk
+            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
+         END DO
+         !
+      ENDIF
 
       ! Local depth and Inverse of the local depth of the water
@@ -475,4 +457,5 @@
       CASE( 'W' )                   !* from T- to W-point : vertical simple mean
          !
+         !zlnwd = 1.0_wp
          pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
          ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
@@ -487,4 +470,5 @@
       CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
          !
+         !zlnwd = 1.0_wp
          pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
@@ -499,4 +483,5 @@
       CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
          !
+         !zlnwd = 1.0_wp
          pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
@@ -676,7 +661,5 @@
                ! Adjust vertical metrics for all wad
                DO jk = 1, jpk
-                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
-                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
+                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )
                END DO
                e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
