Index: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 9342)
+++ branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 9513)
@@ -25,4 +25,8 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE domngb          ! find nearest wet point
+   USE domutl          ! fill closed 3d pool below isf
+   USE domzgr, ONLY : ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat ! import ln_isfsubgl to mask close sea below isf
+   !
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -133,4 +137,5 @@
       INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
       INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
+      INTEGER  ::   jiseed, jjseed           !   -       -
       INTEGER  ::   ios
       INTEGER  ::   isrow                    ! index for ORCA1 starting row
@@ -187,7 +192,4 @@
       END DO  
       
-      ! (ISF) define barotropic mask and mask the ice shelf point
-      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
-      
       DO jk = 1, jpk
          DO jj = 1, jpj
@@ -199,5 +201,33 @@
          END DO  
       END DO  
-
+      !
+      IF ( ln_isfsubgl ) THEN
+         ! check closed wet pool
+         CALL dom_ngb(rn_isfsubgllon, rn_isfsubgllat, jiseed, jjseed, 'T', lwet=.TRUE.)
+         CALL fill_pool( jiseed, jjseed, 1, tmask, -1._wp )
+         ! at this point itab3d (:,1:ijmax,:) can have 3 different values :
+         !              0 where there where already 0
+         !              -1 where the ocean points are connected
+         !              1 where ocean points in tmask are not connected
+         IF (lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*)'dommsk : removal of subglacial lakes '
+            WRITE(numout,*)'~~~~~~~'
+            WRITE(numout,*)'Number of disconected points : ', COUNT(  (tmask(:,:,:) == 1) )
+            WRITE(numout,*)'lon/lat seed to detect main ocean is: ', rn_isfsubgllon, rn_isfsubgllat
+            WRITE(numout,*)'i/j     seed to detect main ocean is: ', jiseed, jjseed
+         END IF
+         DO jk = 1, jpk
+            WHERE (tmask(:,:,jk) > 0 .AND. misfdep(:,:) > 1) tmask(:,:,jk) = 0 ! remove only subglacial lake (ie similar to close sea only below an ice shelf 
+            WHERE (tmask(:,:,jk) < 0) tmask(:,:,jk) = 1 ! restore mask value
+         END DO
+      END IF
+
+      ! (ISF) define barotropic mask and mask the ice shelf point
+      DO jj = 1, jpj
+         DO ji = 1, jpi   ! vector loop
+            ssmask(ji,jj)  = MIN(1._wp,SUM(tmask(ji,jj,:)))
+         END DO
+      END DO
 !!gm  ????
 #if defined key_zdfkpp
Index: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
===================================================================
--- branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90	(revision 9342)
+++ branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90	(revision 9513)
@@ -28,9 +28,9 @@
 CONTAINS
 
-   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )
+   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, lwet, kkk )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dom_ngb  ***
       !!
-      !! ** Purpose :   find the closest grid point from a given lon/lat position
+      !! ** Purpose :   find the closest wet grid point from a given lon/lat position
       !!
       !! ** Method  :   look for minimum distance in cylindrical projection 
@@ -40,9 +40,14 @@
       REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point
       INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point
+      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used
+      LOGICAL         , INTENT(in   ), OPTIONAL :: lwet ! logical to decide if we look for every where (false) 
+                                                                !    or only over the wet point (true)
       CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W'
       !
+      INTEGER :: ik         ! working level
       INTEGER , DIMENSION(2) ::   iloc
       REAL(wp)               ::   zlon, zmini
       REAL(wp), POINTER, DIMENSION(:,:) ::  zglam, zgphi, zmask, zdist
+      LOGICAL :: llwet      ! working logical
       !!--------------------------------------------------------------------
       !
@@ -51,21 +56,46 @@
       CALL wrk_alloc( jpi, jpj, zglam, zgphi, zmask, zdist )
       !
+      ! select lat/lon variable
+      SELECT CASE( cdgrid )
+         CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:)
+         CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:)
+         CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:)
+         CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:)
+      END SELECT
+      !
+      ! select vertical level
+      ik   = 1
+      IF ( PRESENT(kkk) ) ik=kkk
+      !
+      ! select mask variable
       zmask(:,:) = 0._wp
-      SELECT CASE( cdgrid )
-      CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,1)
-      CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,1)
-      CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,1)
-      CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,1)
-      END SELECT
-
-      zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
-      zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
-      IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
-      IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
-
-      zglam(:,:) = zglam(:,:) - zlon
+      llwet = .TRUE.
+      IF ( PRESENT(lwet)) llwet=lwet 
+      IF (llwet) THEN
+         SELECT CASE( cdgrid )
+            CASE( 'U' )  ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik)
+            CASE( 'V' )  ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik)
+            CASE( 'F' )  ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik)
+            CASE DEFAULT ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik)
+         END SELECT
+      ELSE
+         zmask(nldi:nlei,nldj:nlej)=1.e0
+      END IF
+      !
+      ! compute distance
+      IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN
+         zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
+         zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
+         IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
+         IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
+         zglam(:,:) = zglam(:,:) - zlon
+      ELSE
+         zglam(:,:) = zglam(:,:) - plon
+      END IF
+      !
       zgphi(:,:) = zgphi(:,:) - plat
       zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:)
-      
+      !
+      ! find min
       IF( lk_mpp ) THEN  
          CALL mpp_minloc( zdist(:,:), zmask, zmini, kii, kjj)
Index: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domutl.F90
===================================================================
--- branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domutl.F90	(revision 9513)
+++ branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domutl.F90	(revision 9513)
@@ -0,0 +1,268 @@
+MODULE domutl
+   !!
+   !!======================================================================
+   !!                       ***  MODULE  closea  ***
+   !! modutl : specific module to do
+   !!              - 2d flood filling (closea)
+   !!              - 3d flood filling (zgr_isf)
+   !!
+   !!======================================================================
+   !! History :   4.0  !  03-18  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+   !!
+   !!----------------------------------------------------------------------
+   !!   fillpool    : replace all point connected to the seed by a fillvalue 
+   !!----------------------------------------------------------------------
+   !!
+   USE dom_oce         ! ocean space and time domain
+   USE domngb
+   USE in_out_manager    ! I/O manager
+   USE lib_mpp
+   USE lbclnk
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC fill_pool      ! routine called by domain module
+
+   INTERFACE fill_pool
+      MODULE PROCEDURE FillPool2D, FillPool3D
+   END INTERFACE
+
+   TYPE idx
+      INTEGER :: i,j,k
+      TYPE(idx), POINTER :: next
+   END TYPE idx
+
+CONTAINS
+   SUBROUTINE FillPool3D(kiseed, kjseed, kkseed, rdta, rfill)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE FillPool3D  ***
+      !!
+      !! ** Purpose :  Replace all area surrounding by mask value by kifill value
+      !!
+      !! ** Method  :  flood fill algorithm
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER(KIND=4),                  INTENT(in)    :: kiseed, kjseed, kkseed
+      REAL(wp),                         INTENT(in)    :: rfill    ! filling value
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: rdta     ! input data
+      REAL(wp), DIMENSION(jpi,jpj,jpk)                :: rseedmap, rseedmap_b  !
+ 
+      INTEGER :: ip=0                 ! size of the pile
+      INTEGER :: ii  , ij  , ik   , kii, kjj, jj, kk    ! working integer
+      INTEGER :: iip1, ijp1, ikp1
+      INTEGER :: iim1, ijm1, ikm1
+      INTEGER :: nseed
+      TYPE (idx), POINTER :: seed
+      !!---------------------------------------------------------------------- 
+      !
+      ! initialisation seed
+      NULLIFY(seed)
+      !
+      ! create the first seed and update nseed for all processors
+      nseed = 0
+      ii=kiseed; ij=kjseed ; ik=kkseed
+      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
+         ii = mi0(ii) ; ij = mj0(ij);                             ! conversion to local index
+         IF (rdta(ii,ij,ik) > 0 ) THEN                            ! create seed if not on land
+            CALL add_head_idx(seed,ii,ij,ik)
+            rdta    (ii,ij,ik) = rfill
+            rseedmap(ii,ij,ik) = 1.
+         END IF
+      END IF
+      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! nseed =0 means on land => WARNING later on
+      !
+      ! loop until the pile size is 0 or if the pool is larger than the critical size
+      IF (nseed > 0) THEN
+         ! seed on ocean continue
+         DO WHILE ( nseed /= 0 )
+            DO WHILE ( ASSOCIATED(seed) )
+               ip=ip+1
+               ii=seed%i; ij=seed%j ; ik=seed%k ; rseedmap(ii,ij,ik)=1.
+               ! 
+               ! update bathy and update pile size
+               CALL del_idx(seed) 
+               !
+               ! check neighbour cells
+               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
+               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
+               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
+               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
+               ikp1=ik+1 ; IF (ikp1 == jpk+1 ) ikp1=ik
+               ikm1=ik-1 ; IF (ikm1 == 0     ) ikm1=ik
+               !
+               ! ji direction
+               IF (rdta(ii,ijp1,ik) > 0 .AND. rdta(ii,ijp1,ik) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ijp1,ik) ; rdta(ii,ijp1,ik) = rfill ; rseedmap(ii,ijp1,ik)=1. ; END IF
+               IF (rdta(ii,ijm1,ik) > 0 .AND. rdta(ii,ijm1,ik) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ijm1,ik) ; rdta(ii,ijm1,ik) = rfill ; rseedmap(ii,ijm1,ik)=1. ; END IF
+               !
+               ! jj direction
+               IF (rdta(iip1,ij,ik) > 0 .AND. rdta(iip1,ij,ik) /= rfill ) THEN ; CALL add_head_idx(seed,iip1,ij,ik) ; rdta(iip1,ij,ik) = rfill ; rseedmap(iip1,ij,ik)=1. ; END IF
+               IF (rdta(iim1,ij,ik) > 0 .AND. rdta(iim1,ij,ik) /= rfill ) THEN ; CALL add_head_idx(seed,iim1,ij,ik) ; rdta(iim1,ij,ik) = rfill ; rseedmap(iim1,ij,ik)=1. ; END IF
+               !
+               ! jk direction
+               IF (rdta(ii,ij,ikp1) > 0 .AND. rdta(ii,ij,ikp1) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ij,ikp1) ; rdta(ii,ij,ikp1) = rfill ; rseedmap(ii,ij,ik)=1. ; END IF
+               IF (rdta(ii,ij,ikm1) > 0 .AND. rdta(ii,ij,ikm1) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ij,ikm1) ; rdta(ii,ij,ikm1) = rfill ; rseedmap(ii,ij,ik)=1. ; END IF
+            END DO
+            !
+            ! exchange seed
+            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! this is the sum of all the point check this iteration
+            !
+            rseedmap_b(:,:,:)=rseedmap(:,:,:)
+            CALL lbc_lnk(rseedmap,'T',1.)
+            !
+            ! build new list of seed
+            DO ii=1,jpi
+               DO jj=1,jpj
+                  DO kk=1,jpk
+                     IF (rseedmap(ii,jj,kk) > 0.0 .AND. rseedmap(ii,jj,kk) /= rseedmap_b(ii,jj,kk) .AND. rdta(ii,jj,kk) > 0 .AND. rdta(ii,jj,kk) /= rfill) THEN
+                        CALL add_head_idx(seed,ii,jj,kk)
+                     END IF
+                  END DO
+               END DO
+            END DO
+            !
+            ! reset map of seed
+            rseedmap(:,:,:)=0.0
+            !
+         END DO
+      ELSE
+         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do'
+      END IF
+ 
+   END SUBROUTINE FillPool3D
+   SUBROUTINE FillPool2D(kiseed, kjseed, rdta, rfill)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE FillPool3D  ***
+      !!
+      !! ** Purpose :  Replace all area surrounding by mask value by kifill value
+      !!
+      !! ** Method  :  flood fill algorithm
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER(KIND=4),            INTENT(in)    :: kiseed, kjseed
+      REAL(wp),                   INTENT(in)    :: rfill    ! filling value
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: rdta     ! input data
+      REAL(wp), DIMENSION(jpi,jpj) :: rseedmap, rseedmap_b 
+ 
+      INTEGER :: ip=0                     ! size of the pile
+      INTEGER :: ii  , ij  , jj, kii, kjj     ! working integer
+      INTEGER :: iip1, ijp1               ! working integer
+      INTEGER :: iim1, ijm1
+      INTEGER :: nseed
+      TYPE (idx), POINTER :: seed
+      !!----------------------------------------------------------------------
+      !
+      ! initialisation seed
+      NULLIFY(seed)
+      !
+      ! create the first seed and update nseed for all processors
+      nseed = 0
+      ii = kiseed ; ij = kjseed ;
+      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
+         ii = mi0(ii) ; ij = mj0(ij)                              ! conversion to local index
+         IF (rdta(ii,ij) > 0 ) THEN                               ! create seed if not on land
+            CALL add_head_idx(seed,ii,ij,1)
+            rdta    (ii,ij) = rfill
+            rseedmap(ii,ij) = 1.
+         END IF
+      END IF
+      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! nseed =0 means on land => WARNING later on
+      !
+      ! loop until the pile size is 0 or if the pool is larger than the critical size
+      IF (nseed > 0) THEN
+         ! seed on ocean continue
+         DO WHILE ( nseed .NE. 0 )
+            DO WHILE ( ASSOCIATED(seed) )
+               ip=ip+1
+               ii=seed%i; ij=seed%j ; rseedmap(ii,ij)=1.
+               ! update bathy and update pile size
+               CALL del_idx(seed) 
+               ! 
+               ! check neighbour cells
+               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
+               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
+               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
+               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
+               !
+               ! ji direction
+               IF (rdta(ii  ,ijp1) > 0 .AND. rdta(ii  ,ijp1) /= rfill ) THEN ; CALL add_head_idx(seed,ii  ,ijp1,1) ; rdta(ii  ,ijp1) = rfill ; rseedmap(ii  ,ijp1)=1. ; END IF
+               IF (rdta(ii  ,ijm1) > 0 .AND. rdta(ii  ,ijm1) /= rfill ) THEN ; CALL add_head_idx(seed,ii  ,ijm1,1) ; rdta(ii  ,ijm1) = rfill ; rseedmap(ii  ,ijm1)=1. ; END IF
+               !
+               ! jj direction
+               IF (rdta(iip1,ij  ) > 0 .AND. rdta(iip1,ij  ) /= rfill ) THEN ; CALL add_head_idx(seed,iip1,ij  ,1) ; rdta(iip1,ij  ) = rfill ; rseedmap(iip1,ij  )=1. ; END IF
+               IF (rdta(iim1,ij  ) > 0 .AND. rdta(iim1,ij  ) /= rfill ) THEN ; CALL add_head_idx(seed,iim1,ij  ,1) ; rdta(iim1,ij  ) = rfill ; rseedmap(iim1,ij  )=1. ; END IF
+            END DO
+            !
+            ! exchange seed
+            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! this is the sum of all the point check this iteration
+            !
+            rseedmap_b(:,:)=rseedmap(:,:)
+            CALL lbc_lnk(rseedmap,'T',1.)
+            !
+            ! build new list of seed
+            DO ii=1,jpi
+               DO jj=1,jpj
+                  IF (rseedmap(ii,jj) > 0.0 .AND. rseedmap(ii,jj) /= rseedmap_b(ii,jj) .AND. rdta(ii,jj) > 0 .AND. rdta(ii,jj) /= rfill) THEN
+                     CALL add_head_idx(seed, ii, jj, 1)
+                  END IF
+               END DO
+            END DO 
+            !
+            ! reset map of seed
+            rseedmap(:,:)=0.0
+            !
+         END DO
+      ELSE
+         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do'
+      END IF
+ 
+   END SUBROUTINE FillPool2D
+   !
+   !
+   ! subroutine to deals with link list
+   !
+   SUBROUTINE create_idx(pt_idx, ki, kj, kk)
+      TYPE (idx), POINTER :: pt_idx
+      INTEGER, INTENT(in) :: ki, kj, kk
+      !
+      ! initialised all field to NULL()
+      NULLIFY(pt_idx)
+      !
+      ! allocate new element
+      ALLOCATE(pt_idx)
+      pt_idx%i=ki ; pt_idx%j=kj ; pt_idx%k=kk ;
+      pt_idx%next => NULL()  
+   END SUBROUTINE create_idx
+
+   SUBROUTINE add_head_idx(pt_idx, ki, kj, kk)
+      TYPE (idx), POINTER :: pt_idx
+      TYPE (idx), POINTER :: zpt_new
+      INTEGER, INTENT(in) :: ki, kj, kk
+      !
+      ! allocate new element
+      ALLOCATE(zpt_new)
+      zpt_new%i=ki ; zpt_new%j=kj ; zpt_new%k=kk ;
+      zpt_new%next => NULL()
+      !
+      ! linked new element to main list
+      zpt_new%next => pt_idx
+      pt_idx => zpt_new
+   END SUBROUTINE add_head_idx
+
+   SUBROUTINE del_idx(pt_idx)
+      TYPE (idx), POINTER :: pt_idx
+      TYPE (idx), POINTER :: zpt_tmp
+      !
+      ! delete first element
+      IF (ASSOCIATED(pt_idx%next)) THEN
+         zpt_tmp => pt_idx%next
+         DEALLOCATE(pt_idx)
+         NULLIFY(pt_idx)
+         pt_idx => zpt_tmp
+      ELSE
+         NULLIFY(pt_idx)
+      ENDIF
+   END SUBROUTINE del_idx
+   !
+END MODULE domutl
Index: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9342)
+++ branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 9513)
@@ -50,4 +50,7 @@
    PUBLIC   dom_zgr        ! called by dom_init.F90
 
+   ! remove subglacial lake under the ice sheet.
+   LOGICAL , PUBLIC  :: ln_isfsubgl
+   REAL(wp), PUBLIC  :: rn_isfsubgllon, rn_isfsubgllat
    !                              !!* Namelist namzgr_sco *
    LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
@@ -1293,14 +1296,14 @@
       ! NAMELIST
       INTEGER  ::   ios
-      INTEGER  :: nn_kisfmax    = 999. !  
-      REAL(wp) :: rn_isfdep_min = 10.0_wp ! ice shelf minimal thickness 
+      INTEGER  :: nn_kisfmax    = 999.                              !  
+      REAL(wp) :: rn_isfdep_min = 10.0_wp                           ! ice shelf minimal thickness 
       REAL(wp) :: rn_isfhw_deep = 1.e-3 , rn_isfhw_shallow = 1.0e-3 ! threshold to define grounding line in deep/shallow water
-      REAL(wp) :: rn_isfshallow = 0.0_wp  ! threshold to define shallow ice shelf cavity
-      REAL(wp) :: rn_zisfmax    = 6000.0_wp ! maximun meter of ice we are allowed to dig to assure connectivity
-      LOGICAL  :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. !  remove cheminey, assure connectivity
+      REAL(wp) :: rn_isfshallow = 0.0_wp                            ! threshold to define shallow ice shelf cavity
+      REAL(wp) :: rn_zisfmax    = 6000.0_wp                         ! maximun meter of ice we are allowed to dig to assure connectivity
+      LOGICAL  :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. !  remove cheminey, assure connectivity, remove channel in water column (on z space not level space)
       !!---------------------------------------------------------------------
       NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, &
               &           rn_isfdep_min, rn_isfhw_deep, rn_isfhw_shallow, rn_isfshallow, &
-              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel
+              &           ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat
       !
       IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
