Index: NEMO/branches/2021/dev_r14318_RK3_stage1/cfgs/WED025/cpp_WED025.fcm
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/cfgs/WED025/cpp_WED025.fcm	(revision 15155)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/cfgs/WED025/cpp_WED025.fcm	(revision 15157)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys key_xios key_si3
+ bld::tool::fppkeys key_xios key_si3  key_isf
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domzgr.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domzgr.F90	(revision 15155)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domzgr.F90	(revision 15157)
@@ -224,5 +224,5 @@
       INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level
       !
-      INTEGER  ::   jk     ! dummy loop index
+      INTEGER  ::   ji,jj,jk     ! dummy loop index
       INTEGER  ::   inum   ! local logical unit
       REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
@@ -249,4 +249,10 @@
       CALL iom_get( inum, 'ln_isfcav', z_cav )
       IF( z_cav == 0._wp ) THEN   ;   ld_isfcav = .false.   ;   ELSE   ;   ld_isfcav = .true.   ;   ENDIF
+      !
+      !                          !* ocean top and bottom level
+      CALL iom_get( inum, jpdom_global, 'top_level'    , z2d   )   ! 1st wet T-points (ISF)
+      k_top(:,:) = NINT( z2d(:,:) )
+      CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d   )   ! last wet T-points
+      k_bot(:,:) = NINT( z2d(:,:) )
       !
       !                          !* vertical scale factors
@@ -278,4 +284,16 @@
          CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth
          CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths
+#if defined key_qco && key_isf
+!!st         IF( ln_isfcav ) THEN   ! ISF: adjust the calculation of T-point depth done in e3_to_depth subroutine
+!            pdept(:,:,1) = 0.5_wp * pe3w(:,:,1)       ! already computed by e3_to_depth
+            DO_3D( 1, 1, 1, 1, 2, jpk )                     ! vertical sum at partial cell xxxx other level 
+               IF( jk == k_top(ji,jj) ) THEN                    ! first ocean point : partial cell
+                  gdept_0(ji,jj,jk) = gdepw_0(ji,jj,jk  ) + 0.5_wp * e3w_0(ji,jj,jk)   ! = risfdep + 1/2 e3w_0(mikt)
+               ELSE                                            ! other level 
+                  gdept_0(ji,jj,jk) = gdept_0(ji,jj,jk-1) +          e3w_0(ji,jj,jk)
+               ENDIF
+            END_3D
+!!st         ENDIF
+#endif
          IF(lwp) THEN
             WRITE(numout,*)
@@ -285,10 +303,4 @@
          ENDIF
       ENDIF
-      !
-      !                          !* ocean top and bottom level
-      CALL iom_get( inum, jpdom_global, 'top_level'    , z2d   )   ! 1st wet T-points (ISF)
-      k_top(:,:) = NINT( z2d(:,:) )
-      CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d   )   ! last wet T-points
-      k_bot(:,:) = NINT( z2d(:,:) )
       !
       ! reference depth for negative bathy (wetting and drying only)
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynhpg.F90	(revision 15155)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynhpg.F90	(revision 15157)
@@ -554,15 +554,14 @@
       ! compute rhd at the ice/oce interface (ocean side)
       ! usefull to reduce residual current in the test case ISOMIP with no melting
-      DO ji = 1, jpi
-         DO jj = 1, jpj
-            ikt = mikt(ji,jj)
-            zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm)
-            zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm)
-            zdep_top(ji,jj)  = MAX( risfdep(ji,jj) , gdept(ji,jj,1,Kmm) )
-         END DO
-      END DO
-!!st      CALL eos( zts_top, zdep_top, zrhd_top )
-      CALL eos( zts_top, risfdep, zrhd_top )
-
+      DO_2D( 1, 1, 1, 1 )
+         ikt = mikt(ji,jj)
+         zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm)
+         zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm)
+      END_2D
+      DO_2D( 1, 1, 1, 1 )
+         zdep_top(ji,jj)  = MAX( risfdep(ji,jj) , gdept(ji,jj,1,Kmm) )
+      END_2D
+      CALL eos( zts_top, zdep_top, zrhd_top )
+      !
       !                     !===========================!
       !                     !=====  surface value  =====!
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/ISF/isfload.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/ISF/isfload.F90	(revision 15155)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/ISF/isfload.F90	(revision 15157)
@@ -85,6 +85,9 @@
       !
       DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
+#if defined key_qco && key_isf
+         CALL eos( zts_top(:,:,:), gdept_0(:,:,jk), zrhd(:,:,jk) )
+#else
          CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) )
-!!st ==>> CALL eos( zts_top(:,:,:), gdept_0(:,:,jk), zrhd(:,:,jk) )
+#endif
       END DO
       !
@@ -93,10 +96,20 @@
       !
       !                                !- Surface value + ice shelf gradient
-      pload(:,:) = 0._wp                      ! compute pressure due to ice shelf load 
-      DO_2D( 1, 1, 1, 1 )
+      pload(:,:) = 0._wp                      ! compute 2 x pressure due to ice shelf load 
+      DO_2D( 1, 1, 1, 1 )                     ! since is will be multiply by 0.5 in hpg_isf
          ikt = mikt(ji,jj)
          !
          IF ( ikt > 1 ) THEN
             !                                 ! top layer of the ice shelf
+#if defined key_qco && key_isf
+            pload(ji,jj) = pload(ji,jj) + zrhd(ji,jj,1) * e3w_0(ji,jj,1)
+            !
+            DO jk = 2, ikt-1                  ! core layers of the ice shelf
+               pload(ji,jj) = pload(ji,jj) + (zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_0(ji,jj,jk)
+            END DO
+            !                                 ! deepest part of the ice shelf (between deepest T point and ice/ocean interface
+            pload(ji,jj) = pload(ji,jj) + ( zrhdtop_isf(ji,jj) +    zrhd(ji,jj,ikt-1) )   &
+               &                        * (     risfdep(ji,jj) - gdept_0(ji,jj,ikt-1) )
+#else
             pload(ji,jj) = pload(ji,jj)   &
                &         + zrhd (ji,jj,1) * e3w(ji,jj,1,Kmm)
@@ -109,5 +122,5 @@
             pload(ji,jj) = pload(ji,jj) + ( zrhdtop_isf(ji,jj) +  zrhd(ji,jj,ikt-1)     )   &
                &                        * (     risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) )
-!!st ==>>     &                        * (     risfdep(ji,jj) - gdept_0(ji,jj,ikt-1) )
+#endif
             !
          END IF
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traqsr.F90	(revision 15155)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traqsr.F90	(revision 15157)
@@ -786,32 +786,29 @@
       zcoef    =  zprec * rho0_rcp / ( rDt * zQmax * pfr)
       !
-      klev = jpkm1      ! Level of light extinction
-      DO jk = jpkm1, 1, -1
-         IF( SUM( tmask(:,:,jk) ) > 0 ) THEN    ! ocean point at that level
-            zdw  = MAXVAL( gdepw_0(:,:,jk+1) *       wmask(:,:,jk)         )    ! max w-depth at jk+1 level
-            ze3t = MINVAL(   e3t_0(:,:,jk  ) , mask=(wmask(:,:,jk+1)==1)   )    ! minimum e3t at jk   level
-            zhext =  - pL * LOG( zcoef * ze3t )                                 ! extinction depth
-            IF( zdw >= zhext )   klev = jk                                      ! last T-level reached by Qsr
-         ELSE                                   ! only land point at level jk
-            klev = jk                                                           ! local domain sea-bed level 
-         ENDIF
-      END DO
+      IF( ln_zco .OR. ln_zps ) THEN      ! z- or zps coordinate (use 1D ref vertcial coordinate)
+         klev = jpkm1                              ! Level of light extinction zco / zps
+         DO jk = jpkm1, 1, -1
+            zdw  = gdepw_1d(jk+1)                  ! max w-depth at jk+1 level
+            ze3t =   e3t_1d(jk  )                  ! minimum e3t at jk   level
+            zhext =  - pL * LOG( zcoef * ze3t )    ! extinction depth
+            IF( zdw >= zhext )   klev = jk         ! last T-level reached by Qsr
+         END DO
+      ELSE                               ! s- or s-z- coordinate (use 3D vertical coordinate)
+         klev = jpkm1                              ! Level of light extinction 
+         DO jk = jpkm1, 1, -1    ! 
+            IF( SUM( tmask(:,:,jk) ) > 0 ) THEN    ! ocean point at that level
+               zdw  = MAXVAL( gdepw_0(:,:,jk+1) *       wmask(:,:,jk)         )    ! max w-depth at jk+1 level
+               ze3t = MINVAL(   e3t_0(:,:,jk  ) , mask=(wmask(:,:,jk+1)==1)   )    ! minimum e3t at jk   level
+               zhext =  - pL * LOG( zcoef * ze3t )                                 ! extinction depth
+               IF( zdw >= zhext )   klev = jk                                      ! last T-level reached by Qsr
+            ELSE                                   ! only land point at level jk
+               klev = jk                                                           ! local domain sea-bed level 
+            ENDIF
+         END DO
+         CALL mpp_max('tra_qsr', klev)             ! needed for reproducibility   !!st may be modified to avoid this comm.
+         !                                                                        !!st use ssmask to remove the comm ?
+      ENDIF
       !
 !!st      IF(lwp) WRITE(numout,*) '                level of e3t light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m'
-!!st      !
-!!st      klev = jpkm1      ! Level of light extinction
-!!st      DO jk = jpkm1, 1, -1
-!!st         IF( SUM( tmask(:,:,jk) ) > 0 ) THEN    ! ocean point at that level
-!!st            zdw  = MAXVAL( gdepw_0(:,:,jk+1) * wmask(:,:,jk) )    ! max w-depth at jk+1 level
-!!st            ze3t = MINVAL(   e3t_0(:,:,jk  )                 )    ! minimum e3t at jk   level
-!!st            zhext =  - pL * LOG( zcoef * pL )                   ! extinction depth
-!!st            IF( zdw >= zhext )   klev = jk                        ! last T-level reached by Qsr
-!!st         ELSE                                   ! only land point at level jk
-!!st            klev = jk                                             ! local domain sea-bed level 
-!!st         ENDIF
-!!st      END DO
-!!st      !
-!!st      IF(lwp) WRITE(numout,*) '                level of pL light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m'
-      !
    END FUNCTION qsr_ext_lev
 
Index: NEMO/branches/2021/dev_r14318_RK3_stage1/tests/ISOMIP+/cpp_ISOMIP+.fcm
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/tests/ISOMIP+/cpp_ISOMIP+.fcm	(revision 15155)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/tests/ISOMIP+/cpp_ISOMIP+.fcm	(revision 15157)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys   key_xios 
+ bld::tool::fppkeys   key_xios key_qco key_isf 
