Index: branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5473)
+++ branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5575)
@@ -17,4 +17,5 @@
    !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
    !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case  
+   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye  
    !!----------------------------------------------------------------------
 
@@ -35,5 +36,4 @@
    USE oce               ! ocean variables
    USE dom_oce           ! ocean domain
-   USE sbc_oce           ! surface variable (isf)
    USE closea            ! closed seas
    USE c1d               ! 1D vertical configuration
@@ -298,15 +298,17 @@
       ENDIF
 
+      IF ( ln_isfcav ) THEN
 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
-      DO jk = 1, jpkm1
-         e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
-      END DO
-      e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
-
-      DO jk = 2, jpk
-         e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
-      END DO
-      e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
+         DO jk = 1, jpkm1
+            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
+         END DO
+         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
+
+         DO jk = 2, jpk
+            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
+         END DO
+         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
+      END IF
 
 !!gm BUG in s-coordinate this does not work!
@@ -471,26 +473,4 @@
          misfdep(:,:)=1
          !
-         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
-         IF( cp_cfg == "isomip" ) THEN 
-	        ! 
-	        risfdep(:,:)=200.e0 
-	        misfdep(:,:)=1 
-	        ij0 = 1 ; ij1 = 40 
-	        DO jj = mj0(ij0), mj1(ij1) 
-	           risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
-                END DO 
-            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
-	        ! 
-         ELSEIF ( cp_cfg == "isomip2" ) THEN
-         ! 
-            risfdep(:,:)=0.e0
-            misfdep(:,:)=1
-            ij0 = 1 ; ij1 = 40
-            DO jj = mj0(ij0), mj1(ij1)
-               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
-            END DO
-            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
-         END IF
-         !
          DEALLOCATE( idta, zdta )
          !
@@ -534,7 +514,11 @@
          IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
             CALL iom_open ( 'bathy_meter.nc', inum ) 
-            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
+            IF ( ln_isfcav ) THEN
+               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
+            ELSE
+               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  )
+            END IF
             CALL iom_close( inum )
-            !  
+            !                                                
             risfdep(:,:)=0._wp         
             misfdep(:,:)=1             
@@ -584,7 +568,9 @@
       IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
          ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
-         WHERE (bathy == risfdep)
-            bathy   = 0.0_wp ; risfdep = 0.0_wp
-         END WHERE
+         IF ( ln_isfcav ) THEN
+            WHERE (bathy == risfdep)
+               bathy   = 0.0_wp ; risfdep = 0.0_wp
+            END WHERE
+         END IF
          ! end patch
          IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
@@ -961,4 +947,312 @@
       !!----------------------------------------------------------------------
       !!
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
+      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
+      LOGICAL  ::   ll_print         ! Allow  control print for debugging
+      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
+      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
+      REAL(wp) ::   zmax             ! Maximum depth
+      REAL(wp) ::   zdiff            ! temporary scalar
+      REAL(wp) ::   zrefdep          ! temporary scalar
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zprt )
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
+      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
+      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
+
+      ll_print = .FALSE.                   ! Local variable for debugging
+      
+      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
+         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
+      ENDIF
+
+
+      ! bathymetry in level (from bathy_meter)
+      ! ===================
+      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
+      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
+      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
+      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
+      END WHERE
+
+      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
+      ! find the number of ocean levels such that the last level thickness
+      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
+      ! e3t_1d is the reference level thickness
+      DO jk = jpkm1, 1, -1
+         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
+      END DO
+
+      IF ( ln_isfcav ) CALL zgr_isf
+
+      ! Scale factors and depth at T- and W-points
+      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
+         gdept_0(:,:,jk) = gdept_1d(jk)
+         gdepw_0(:,:,jk) = gdepw_1d(jk)
+         e3t_0  (:,:,jk) = e3t_1d  (jk)
+         e3w_0  (:,:,jk) = e3w_1d  (jk)
+      END DO
+      ! 
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ik = mbathy(ji,jj)
+            IF( ik > 0 ) THEN               ! ocean point only
+               ! max ocean level case
+               IF( ik == jpkm1 ) THEN
+                  zdepwp = bathy(ji,jj)
+                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
+                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
+                  e3t_0(ji,jj,ik  ) = ze3tp
+                  e3t_0(ji,jj,ik+1) = ze3tp
+                  e3w_0(ji,jj,ik  ) = ze3wp
+                  e3w_0(ji,jj,ik+1) = ze3tp
+                  gdepw_0(ji,jj,ik+1) = zdepwp
+                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
+                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
+                  !
+               ELSE                         ! standard case
+                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
+                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
+                  ENDIF
+!gm Bug?  check the gdepw_1d
+                  !       ... on ik
+                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
+                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
+                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
+                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
+                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
+                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
+                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
+                  !       ... on ik+1
+                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
+               ENDIF
+            ENDIF
+         END DO
+      END DO
+      !
+      it = 0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ik = mbathy(ji,jj)
+            IF( ik > 0 ) THEN               ! ocean point only
+               e3tp (ji,jj) = e3t_0(ji,jj,ik)
+               e3wp (ji,jj) = e3w_0(ji,jj,ik)
+               ! test
+               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
+               IF( zdiff <= 0._wp .AND. lwp ) THEN 
+                  it = it + 1
+                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
+                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
+                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
+                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
+               ENDIF
+            ENDIF
+         END DO
+      END DO
+      !
+      IF ( ln_isfcav ) THEN
+      ! (ISF) Definition of e3t, u, v, w for ISF case
+         DO jj = 1, jpj 
+            DO ji = 1, jpi 
+               ik = misfdep(ji,jj) 
+               IF( ik > 1 ) THEN               ! ice shelf point only 
+                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
+                  gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
+!gm Bug?  check the gdepw_0 
+               !       ... on ik 
+                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
+                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
+                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
+                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
+                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
+
+                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
+                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
+                  ENDIF 
+               !       ... on ik / ik-1 
+                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
+                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
+! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
+                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
+               ENDIF 
+            END DO 
+         END DO 
+      ! 
+         it = 0 
+         DO jj = 1, jpj 
+            DO ji = 1, jpi 
+               ik = misfdep(ji,jj) 
+               IF( ik > 1 ) THEN               ! ice shelf point only 
+                  e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
+                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
+               ! test 
+                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
+                  IF( zdiff <= 0. .AND. lwp ) THEN  
+                     it = it + 1 
+                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
+                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
+                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
+                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
+                  ENDIF 
+               ENDIF 
+            END DO 
+         END DO 
+      END IF
+      ! END (ISF)
+
+      ! Scale factors and depth at U-, V-, UW and VW-points
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
+         e3u_0 (:,:,jk) = e3t_1d(jk)
+         e3v_0 (:,:,jk) = e3t_1d(jk)
+         e3uw_0(:,:,jk) = e3w_1d(jk)
+         e3vw_0(:,:,jk) = e3w_1d(jk)
+      END DO
+      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
+               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
+               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
+               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      IF ( ln_isfcav ) THEN
+      ! (ISF) define e3uw (adapted for 2 cells in the water column)
+         DO jj = 2, jpjm1 
+            DO ji = 2, fs_jpim1   ! vector opt. 
+               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
+               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
+               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
+               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
+            END DO
+         END DO
+      END IF
+
+      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
+      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
+      !
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
+         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
+      END DO
+      
+      ! Scale factor at F-point
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
+         e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
+      !
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+!!gm  bug ? :  must be a do loop with mj0,mj1
+      ! 
+      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
+      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
+      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
+      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
+      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
+
+      ! Control of the sign
+      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
+      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
+      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
+      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
+     
+      ! Compute gdep3w_0 (vertical sum of e3w)
+      IF ( ln_isfcav ) THEN ! if cavity
+         WHERE (misfdep == 0) misfdep = 1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
+               DO jk = 2, misfdep(ji,jj)
+                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+               END DO
+               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
+               DO jk = misfdep(ji,jj) + 1, jpk
+                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+      ELSE ! no cavity
+         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
+         DO jk = 2, jpk
+            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
+         END DO
+      END IF
+      !                                               ! ================= !
+      IF(lwp .AND. ll_print) THEN                     !   Control print   !
+         !                                            ! ================= !
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik = MAX( mbathy(ji,jj), 1 )
+               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
+               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
+               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
+               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
+               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
+               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
+            END DO
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+      ENDIF  
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
+      !
+   END SUBROUTINE zgr_zps
+
+   SUBROUTINE zgr_isf
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_isf  ***
+      !!   
+      !! ** Purpose :   check the bathymetry in levels
+      !!   
+      !! ** Method  :   THe water column have to contained at least 2 cells
+      !!                Bathymetry and isfdraft are modified (dig/close) to respect
+      !!                this criterion.
+      !!                 
+      !!   
+      !! ** Action  : - test compatibility between isfdraft and bathy 
+      !!              - bathy and isfdraft are modified
+      !!----------------------------------------------------------------------
+      !!   
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       INTEGER  ::   ik, it           ! temporary integers
@@ -971,45 +1265,15 @@
       REAL(wp) ::   zdiff            ! temporary scalar
       REAL(wp) ::   zrefdep          ! temporary scalar
-      REAL(wp) ::   zbathydiff, zrisfdepdiff 
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH)
-      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
+      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
+      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('zgr_zps')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zprt )
+      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
+      !
       CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
-      CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep)
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps'
-      IF(lwp) WRITE(numout,*) '    ~~~~~~~ '
-      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
-
-      ll_print = .FALSE.                   ! Local variable for debugging
-      
-      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
-         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
-      ENDIF
-
-      ! bathymetry in level (from bathy_meter)
-      ! ===================
-      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
-      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
-      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
-      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
-      END WHERE
-
-      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
-      ! find the number of ocean levels such that the last level thickness
-      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
-      ! e3t_1d is the reference level thickness
-      DO jk = jpkm1, 1, -1
-         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
-         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
-      END DO
+      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
+
+
       ! (ISF) compute misfdep
       WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1  
@@ -1055,10 +1319,10 @@
             misfdep(jpi,:) = misfdep(  2  ,:) 
          ENDIF
- 
+
          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
             mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
             mbathy(jpi,:) = mbathy(  2  ,:)
          ENDIF
- 
+
          ! split last cell if possible (only where water column is 2 cell or less)
          DO jk = jpkm1, 1, -1
@@ -1078,5 +1342,5 @@
             END WHERE
          END DO
- 
+
  
  ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
@@ -1254,5 +1518,5 @@
  
  ! remove single point "bay" on isf coast line in the ice shelf draft'
-         DO jk = 1, jpk
+         DO jk = 2, jpk
             WHERE (misfdep==0) misfdep=jpk
             zmask=0
@@ -1359,5 +1623,5 @@
                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == 0 ) THEN
+               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
                END IF
@@ -1475,238 +1739,10 @@
       ENDIF 
 
-      ! Scale factors and depth at T- and W-points
-      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
-         gdept_0(:,:,jk) = gdept_1d(jk)
-         gdepw_0(:,:,jk) = gdepw_1d(jk)
-         e3t_0  (:,:,jk) = e3t_1d  (jk)
-         e3w_0  (:,:,jk) = e3w_1d  (jk)
-      END DO
-      ! 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ik = mbathy(ji,jj)
-            IF( ik > 0 ) THEN               ! ocean point only
-               ! max ocean level case
-               IF( ik == jpkm1 ) THEN
-                  zdepwp = bathy(ji,jj)
-                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
-                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
-                  e3t_0(ji,jj,ik  ) = ze3tp
-                  e3t_0(ji,jj,ik+1) = ze3tp
-                  e3w_0(ji,jj,ik  ) = ze3wp
-                  e3w_0(ji,jj,ik+1) = ze3tp
-                  gdepw_0(ji,jj,ik+1) = zdepwp
-                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
-                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
-                  !
-               ELSE                         ! standard case
-                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
-                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
-                  ENDIF
-!gm Bug?  check the gdepw_1d
-                  !       ... on ik
-                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
-                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
-                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
-                  e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
-                     &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
-                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
-                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
-                  !       ... on ik+1
-                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
-                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
-                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
-               ENDIF
-            ENDIF
-         END DO
-      END DO
-      !
-      it = 0
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ik = mbathy(ji,jj)
-            IF( ik > 0 ) THEN               ! ocean point only
-               e3tp (ji,jj) = e3t_0(ji,jj,ik)
-               e3wp (ji,jj) = e3w_0(ji,jj,ik)
-               ! test
-               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
-               IF( zdiff <= 0._wp .AND. lwp ) THEN 
-                  it = it + 1
-                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
-                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
-                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
-                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
-               ENDIF
-            ENDIF
-         END DO
-      END DO
-      !
-      ! (ISF) Definition of e3t, u, v, w for ISF case
-      DO jj = 1, jpj 
-         DO ji = 1, jpi 
-            ik = misfdep(ji,jj) 
-            IF( ik > 1 ) THEN               ! ice shelf point only 
-               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
-               gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
-!gm Bug?  check the gdepw_0 
-               !       ... on ik 
-               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
-                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
-                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
-               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
-               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
-
-               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
-                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
-                ENDIF 
-               !       ... on ik / ik-1 
-               e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
-               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
-! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
-               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
-            ENDIF 
-         END DO 
-      END DO 
-      ! 
-      it = 0 
-      DO jj = 1, jpj 
-         DO ji = 1, jpi 
-            ik = misfdep(ji,jj) 
-            IF( ik > 1 ) THEN               ! ice shelf point only 
-               e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
-               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
-               ! test 
-               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
-               IF( zdiff <= 0. .AND. lwp ) THEN  
-                  it = it + 1 
-                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
-                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
-                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
-                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
-               ENDIF 
-            ENDIF 
-         END DO 
-      END DO 
-      ! END (ISF)
-
-      ! Scale factors and depth at U-, V-, UW and VW-points
-      DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3u_0 (:,:,jk) = e3t_1d(jk)
-         e3v_0 (:,:,jk) = e3t_1d(jk)
-         e3uw_0(:,:,jk) = e3w_1d(jk)
-         e3vw_0(:,:,jk) = e3w_1d(jk)
-      END DO
-      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
-               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
-               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
-               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
-            END DO
-         END DO
-      END DO
-      ! (ISF) define e3uw
-      DO jk = 2,jpk                          
-         DO jj = 1, jpjm1 
-            DO ji = 1, fs_jpim1   ! vector opt. 
-               e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) &
-                 &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) )
-               e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) &
-                 &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) )
-            END DO 
-         END DO 
-      END DO
-      !End (ISF)
-      
-      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
-      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
-      !
-      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
-         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
-         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
-         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
-      END DO
-      
-      ! Scale factor at F-point
-      DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3f_0(:,:,jk) = e3t_1d(jk)
-      END DO
-      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
-      !
-      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
-      END DO
-!!gm  bug ? :  must be a do loop with mj0,mj1
-      ! 
-      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
-      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
-      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
-      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
-      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
-
-      ! Control of the sign
-      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
-      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
-      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
-      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
-     
-      ! Compute gdep3w_0 (vertical sum of e3w)
-      WHERE (misfdep == 0) misfdep = 1
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
-            DO jk = 2, misfdep(ji,jj)
-               gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
-            END DO
-            IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
-            DO jk = misfdep(ji,jj) + 1, jpk
-               gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
-            END DO
-        END DO
-      END DO
-      !                                               ! ================= !
-      IF(lwp .AND. ll_print) THEN                     !   Control print   !
-         !                                            ! ================= !
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX( mbathy(ji,jj), 1 )
-               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
-               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
-               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
-               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
-               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
-               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-      ENDIF  
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
       CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
       CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
-      !
-   END SUBROUTINE zgr_zps
+
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
+      
+   END SUBROUTINE
 
    SUBROUTINE zgr_sco
