Index: utils/tools/NESTING/src/agrif_readwrite.f90
===================================================================
--- utils/tools/NESTING/src/agrif_readwrite.f90	(revision 10160)
+++ utils/tools/NESTING/src/agrif_readwrite.f90	(revision 10381)
@@ -209,5 +209,5 @@
     TYPE(Coordinates) :: Grid
     !
-    CALL read_ncdf_var('mbathy',name,Grid%Bathy_level)    
+    CALL read_ncdf_var(parent_level_name,name,Grid%Bathy_level)    
     !
     WRITE(*,*) ' '
@@ -239,11 +239,11 @@
     CALL write_ncdf_dim(dimnames(2),name,nyfin)
     !      
-    CALL write_ncdf_var('nav_lon',dimnames,name,Grid%nav_lon    ,'float')
-    CALL write_ncdf_var('nav_lat',dimnames,name,Grid%nav_lat    ,'float')
-    CALL write_ncdf_var('mbathy' ,dimnames,name,Grid%bathy_level,'float')
-    !
-    CALL copy_ncdf_att('nav_lon',TRIM(parent_bathy_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
-    CALL copy_ncdf_att('nav_lat',TRIM(parent_bathy_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
-    CALL copy_ncdf_att('mbathy' ,TRIM(parent_bathy_file),name)       
+    CALL write_ncdf_var('nav_lon'        ,dimnames,name,Grid%nav_lon    ,'float')
+    CALL write_ncdf_var('nav_lat'        ,dimnames,name,Grid%nav_lat    ,'float')
+    CALL write_ncdf_var(parent_level_name,dimnames,name,Grid%bathy_level,'float')
+    !
+    CALL copy_ncdf_att('nav_lon'        ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
+    CALL copy_ncdf_att('nav_lat'        ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
+    CALL copy_ncdf_att(parent_level_name,TRIM(parent_bathy_level),name)       
     !
     WRITE(*,*) ' '
@@ -256,8 +256,8 @@
   !
   !*****************************************************
-  !   function Read_Bathy_meter(name,CoarseGrid,ChildGrid)
-  !*****************************************************
-  !
-  INTEGER FUNCTION Read_Bathy_meter(name,CoarseGrid,ChildGrid,Pacifique)
+  !   function read_bathy_coord(name,CoarseGrid,ChildGrid)
+  !*****************************************************
+  !
+  INTEGER FUNCTION read_bathy_coord(name,CoarseGrid,ChildGrid,Pacifique)
     !
     USE io_netcdf
@@ -284,5 +284,5 @@
        CALL read_ncdf_var('nav_lon',name,CoarseGrid%nav_lon)
        CALL read_ncdf_var('nav_lat',name,CoarseGrid%nav_lat)
-       CALL read_ncdf_var(parent_batmet_name,name,CoarseGrid%Bathy_meter)
+       CALL read_ncdf_var(parent_meter_name,name,CoarseGrid%Bathy_meter)
        !            
        IF ( PRESENT(Pacifique) ) THEN
@@ -294,5 +294,5 @@
        ENDIF
        !      
-       Read_Bathy_meter = 1
+       read_bathy_coord = 1
        RETURN      
     ELSE
@@ -411,15 +411,15 @@
     CoarseGrid%Bathy_meter(:,:) = -1.0 * CoarseGrid%Bathy_meter(:,:)
     !      
-    Read_Bathy_meter = 1
+    read_bathy_coord = 1
     RETURN
     !      
-  END FUNCTION Read_Bathy_meter
+  END FUNCTION read_bathy_coord
   ! 
   !
   !*****************************************************
-  !   function Read_Bathy_meter(name,CoarseGrid,ChildGrid)
-  !*****************************************************
-  !
-  INTEGER FUNCTION Read_Bathymeter(name,Grid)
+  !   function read_bathy_meter(name,CoarseGrid,ChildGrid)
+  !*****************************************************
+  !
+  INTEGER FUNCTION read_bathy_meter(name,Grid)
     !
     !
@@ -429,5 +429,5 @@
     TYPE(Coordinates) :: Grid
     !
-    CALL read_ncdf_var(parent_batmet_name,name,Grid%Bathy_meter)    
+    CALL read_ncdf_var(parent_meter_name,name,Grid%Bathy_meter)    
     !
     WRITE(*,*) ' '
@@ -435,7 +435,7 @@
     WRITE(*,*) ' '      
     !
-    Read_Bathymeter = 1
-    !      
-  END FUNCTION Read_Bathymeter
+    read_bathy_meter = 1
+    !      
+  END FUNCTION read_bathy_meter
   !     
   !*****************************************************
@@ -463,12 +463,12 @@
     CALL write_ncdf_dim(dimnames(2),name,ny)
     !     
-    CALL write_ncdf_var('nav_lon'         ,dimnames,name,Grid%nav_lon    ,'float')
-    CALL write_ncdf_var('nav_lat'         ,dimnames,name,Grid%nav_lat    ,'float')
-    CALL write_ncdf_var(parent_batmet_name,dimnames,name,Grid%bathy_meter,'float')
-    CALL write_ncdf_var('weight'          ,dimnames,name,Grid%wgt        ,'float')
-    !
-    CALL copy_ncdf_att('nav_lon'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
-    CALL copy_ncdf_att('nav_lat'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))  
-    CALL copy_ncdf_att(parent_batmet_name,TRIM(parent_bathy_meter),name)
+    CALL write_ncdf_var('nav_lon'        ,dimnames,name,Grid%nav_lon    ,'float')
+    CALL write_ncdf_var('nav_lat'        ,dimnames,name,Grid%nav_lat    ,'float')
+    CALL write_ncdf_var(parent_meter_name,dimnames,name,Grid%bathy_meter,'float')
+    CALL write_ncdf_var('weight'         ,dimnames,name,Grid%wgt        ,'float')
+    !
+    CALL copy_ncdf_att('nav_lon'        ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
+    CALL copy_ncdf_att('nav_lat'        ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))  
+    CALL copy_ncdf_att(parent_meter_name,TRIM(parent_bathy_meter),name)
     !
     WRITE(*,*) ' '
@@ -500,7 +500,7 @@
      !
      CHARACTER(len=1),     DIMENSION(3)     ::   dimnames
-     REAL*8 ,              DIMENSION(N)     ::   e3t_1d, e3w_1d, gdept_1d
-     REAL*8 , ALLOCATABLE, DIMENSION(:,:)   ::   ff_t, ff_f
-     INTEGER, ALLOCATABLE, DIMENSION(:,:)   ::   bottom_level, top_level, zbathy
+     REAL*8 ,              DIMENSION(N)     ::   e3t_1d, e3w_1d, gdept_1d, gdepw_1d
+     REAL*8 , ALLOCATABLE, DIMENSION(:,:)   ::   ff_t, ff_f, zbathy
+     INTEGER, ALLOCATABLE, DIMENSION(:,:)   ::   bottom_level, top_level, mbathy
      REAL*8 , ALLOCATABLE, DIMENSION(:,:,:) ::   e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0
      !
@@ -510,6 +510,6 @@
 
      ! allocate needed arrays for domain_cfg
-     ALLOCATE( ff_t(nx,ny), ff_f(nx,ny) )
-     ALLOCATE( bottom_level(nx,ny), top_level(nx,ny), zbathy(nx,ny) )
+     ALLOCATE( ff_t(nx,ny), ff_f(nx,ny), zbathy(nx,ny) )
+     ALLOCATE( bottom_level(nx,ny), top_level(nx,ny), mbathy(nx,ny) )
      ALLOCATE( e3t_0(nx,ny,N), e3u_0 (nx,ny,N), e3v_0 (nx,ny,N), e3f_0(nx,ny,N),  &
         &      e3w_0(nx,ny,N), e3uw_0(nx,ny,N), e3vw_0(nx,ny,N) )
@@ -527,19 +527,18 @@
      ff_t(:,:) = 2. * omega * SIN( rad * Grid%gphit(:,:) )     !    -        -       -    at t-point
 
-     ! top/bottom levels
-     bottom_level(:,:) = Grid%bathy_level(:,:)
-     top_level(:,:)    = MIN( 1, bottom_level(:,:) )
+     ! bathy in meters
+     zbathy(:,:) = Grid%bathy_meter
 
      ! vertical scale factors
-     CALL zgr_z( e3t_1d, e3w_1d, gdept_1d )    
-     DO jk = 1, N
-        e3t_0  (:,:,jk) = e3t_1d  (jk)
-        e3u_0  (:,:,jk) = e3t_1d  (jk)
-        e3v_0  (:,:,jk) = e3t_1d  (jk)
-        e3f_0  (:,:,jk) = e3t_1d  (jk)
-        e3w_0  (:,:,jk) = e3w_1d  (jk)
-        e3uw_0 (:,:,jk) = e3w_1d  (jk)
-        e3vw_0 (:,:,jk) = e3w_1d  (jk)
-     END DO
+     CALL zgr_z( e3t_1d, e3w_1d, gdept_1d, gdepw_1d )    
+!     DO jk = 1, N
+!        e3t_0  (:,:,jk) = e3t_1d  (jk)
+!        e3u_0  (:,:,jk) = e3t_1d  (jk)
+!        e3v_0  (:,:,jk) = e3t_1d  (jk)
+!        e3f_0  (:,:,jk) = e3t_1d  (jk)
+!        e3w_0  (:,:,jk) = e3w_1d  (jk)
+!        e3uw_0 (:,:,jk) = e3w_1d  (jk)
+!        e3vw_0 (:,:,jk) = e3w_1d  (jk)
+!     END DO
 
      ! logicals and others
@@ -549,4 +548,6 @@
         ln_zps = 1
         ln_zco = 0
+        CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &        ! input
+        &             mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )  ! output
      ELSE
         ln_zps = 0
@@ -554,4 +555,8 @@
      ENDIF
 
+     ! top/bottom levels
+     bottom_level(:,:) = mbathy(:,:)
+     top_level(:,:)    = MIN( 1, mbathy(:,:) )
+
      ! closed domain (agrif)
      jperio = 0
@@ -567,5 +572,4 @@
      top_level(  nx,1:ny) = 0
 
-     zbathy(:,:) = Grid%bathy_meter
      zbathy(1:nx,1   ) = 0.
      zbathy(1:nx,  ny) = 0.
@@ -576,5 +580,5 @@
      ! write domain_cfg.nc
      ! -------------------
-     status = nf90_create(name,NF90_WRITE,ncid)
+     status = nf90_create(name,IOR(NF90_64BIT_OFFSET,NF90_WRITE),ncid)
      status = nf90_close(ncid)
      !
@@ -664,6 +668,6 @@
      WRITE(*,*) ' '
      !
-     DEALLOCATE( ff_t, ff_f )
-     DEALLOCATE( bottom_level, top_level, zbathy )
+     DEALLOCATE( ff_t, ff_f, zbathy )
+     DEALLOCATE( bottom_level, top_level, mbathy )
      DEALLOCATE( e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )
      !      
@@ -672,5 +676,5 @@
   END FUNCTION write_domcfg
   !
-  SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d)
+  SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d, gdepw_1d )
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE zgr_z (from NEMO4) ***
@@ -694,5 +698,4 @@
      !!----------------------------------------------------------------------
      INTEGER  ::   jk                   ! dummy loop indices
-     INTEGER  ::   jpk
      REAL*8 ::   zt, zw                 ! temporary scalars
      REAL*8 ::   zsur, za0, za1, zkth   ! Values set from parameters in
@@ -700,6 +703,5 @@
      REAL*8 ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters 
 
-     REAL*8, DIMENSION(:), INTENT(inout) ::  e3t_1d, e3w_1d, gdept_1d
-     REAL*8, DIMENSION(N)                ::  gdepw_1d
+     REAL*8, DIMENSION(:), INTENT(out) ::  e3t_1d, e3w_1d, gdept_1d, gdepw_1d
      !!----------------------------------------------------------------------
      !
@@ -788,4 +790,240 @@
   END SUBROUTINE zgr_z
 
+  SUBROUTINE zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &
+     &                mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zgr_zps (from NEMO4)  ***
+      !!                     
+      !! ** Purpose :   the depth and vertical scale factor in partial step
+      !!              reference z-coordinate case
+      !!
+      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
+      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
+      !!      a partial step representation of bottom topography.
+      !!
+      !!        The reference depth of model levels is defined from an analytical
+      !!      function the derivative of which gives the reference vertical
+      !!      scale factors.
+      !!        From  depth and scale factors reference, we compute there new value
+      !!      with partial steps  on 3d arrays ( i, j, k ).
+      !!
+      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
+      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
+      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
+      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
+      !!
+      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
+      !!      we find the mbathy index of the depth at each grid point.
+      !!      This leads us to three cases:
+      !!
+      !!              - bathy = 0 => mbathy = 0
+      !!              - 1 < mbathy < jpkm1    
+      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1  
+      !!
+      !!        Then, for each case, we find the new depth at t- and w- levels
+      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw- 
+      !!      and f-points.
+      !! 
+      !!        This routine is given as an example, it must be modified
+      !!      following the user s desiderata. nevertheless, the output as
+      !!      well as the way to compute the model levels and scale factors
+      !!      must be respected in order to insure second order accuracy
+      !!      schemes.
+      !!
+      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
+      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
+      !!      
+      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
+      INTEGER  ::   ik
+      INTEGER  ::   ikb, ikt       ! temporary integers
+      REAL*8 ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
+      REAL*8 ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
+      REAL*8 ::   zmax             ! temporary scalar
+      !
+      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) ::  gdept_0, gdepw_0, zprt
+      !
+      INTEGER,                   INTENT(in   ) ::  nx, ny
+      REAL*8 , DIMENSION(:)    , INTENT(in   ) ::  gdept_1d, gdepw_1d, e3t_1d, e3w_1d
+      REAL*8 , DIMENSION(:,:)  , INTENT(inout) ::  zbathy
+      INTEGER, DIMENSION(:,:)  , INTENT(  out) ::  mbathy
+      REAL*8 , DIMENSION(:,:,:), INTENT(  out) ::  e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0
+      !
+      !!---------------------------------------------------------------------
+      
+      ALLOCATE( zprt(nx,ny,N), gdept_0(nx,ny,N), gdepw_0(nx,ny,N) )
+      !
+      ! bathymetry in level (from bathy_meter)
+      ! ===================
+      zmax = gdepw_1d(N) + e3t_1d(N)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) )
+      zbathy(:,:) = MIN( zmax ,  zbathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
+      WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
+      ELSEWHERE                    ;   mbathy(:,:) = N-1   ! 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 = N-1, 1, -1
+         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+         WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
+      END DO
+
+      ! Scale factors and depth at T- and W-points
+      DO jk = 1, N                        ! 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
+      
+      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
+      !!clem: not implemented yet
+      !!      IF ( ln_isfcav == 1 ) CALL zgr_isf
+      !
+      ! Scale factors and depth at T- and W-points
+!      IF ( ln_isfcav == 0 ) THEN
+         DO jj = 1, ny
+            DO ji = 1, nx
+               ik = mbathy(ji,jj)
+               IF( ik > 0 ) THEN               ! ocean point only
+                  ! max ocean level case
+                  IF( ik == N-1 ) THEN
+                     zdepwp = zbathy(ji,jj)
+                     ze3tp  = zbathy(ji,jj) - gdepw_1d(ik)
+                     ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( 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( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = zbathy(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 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * 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
+         !
+!         DO jj = 1, ny
+!            DO ji = 1, nx
+!               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)
+!               ENDIF
+!            END DO
+!         END DO
+!      END IF
+      !
+      ! Scale factors and depth at U-, V-, UW and VW-points
+      DO jk = 1, N                        ! 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,N                         ! Computed as the minimum of neighbooring scale factors
+         DO jj = 1, ny-1
+            DO ji = 1, nx-1   ! 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 == 1 ) THEN
+!      ! (ISF) define e3uw (adapted for 2 cells in the water column)
+!         DO jj = 2, ny-1 
+!            DO ji = 2, nx-1   ! 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
+      !
+
+      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3u_0 (:,:,jk) == 0. )   e3u_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3v_0 (:,:,jk) == 0. )   e3v_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3uw_0(:,:,jk) == 0. )   e3uw_0(:,:,jk) = e3w_1d(jk)
+         WHERE( e3vw_0(:,:,jk) == 0. )   e3vw_0(:,:,jk) = e3w_1d(jk)
+      END DO
+      
+      ! Scale factor at F-point
+      DO jk = 1, N                        ! initialisation to z-scale factors
+         e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+      DO jk = 1, N                        ! Computed as the minimum of neighbooring V-scale factors
+         DO jj = 1, ny-1
+            DO ji = 1, nx-1   ! 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
+      !
+      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3f_0(:,:,jk) == 0. )   e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+
+      ! Control of the sign
+      IF( MINVAL( e3t_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3t_0 <= 0' 
+      IF( MINVAL( e3w_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3w_0 <= 0' 
+      IF( MINVAL( gdept_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdept_0 <  0' 
+      IF( MINVAL( gdepw_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdepw_0 <  0' 
+     
+      ! Compute gde3w_0 (vertical sum of e3w)
+!      IF ( ln_isfcav ==1 ) THEN ! if cavity
+!         WHERE( misfdep == 0 )   misfdep = 1
+!         DO jj = 1,ny
+!            DO ji = 1,nx
+!               gde3w_0(ji,jj,1) = 0.5 * e3w_0(ji,jj,1)
+!               DO jk = 2, misfdep(ji,jj)
+!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+!               END DO
+!               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5 * e3w_0(ji,jj,misfdep(ji,jj))
+!               DO jk = misfdep(ji,jj) + 1, N
+!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+!               END DO
+!            END DO
+!         END DO
+!      ELSE ! no cavity
+!         gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
+!         DO jk = 2, N
+!            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
+!         END DO
+!      END IF
+      !
+      DEALLOCATE( zprt, gdept_0, gdepw_0 )
+      !
+   END SUBROUTINE zgr_zps
 
   !*****************************************************
@@ -828,5 +1066,5 @@
     LOGICAL :: find        
     i=1
-    DO WHILE ( TRIM(VAR_INTERP(i)) .NE. 'NULL' )     
+    DO WHILE ( TRIM(VAR_INTERP(i)) /= '' )     
        pos = INDEX( TRIM(VAR_INTERP(i)) , TRIM(varname) )
        IF ( pos .NE. 0 ) THEN      
