Index: branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4244)
+++ branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4245)
@@ -374,7 +374,21 @@
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
-            idta(:,:) = jpkm1                            ! before last level
-            zdta(:,:) = gdepw_0(jpk)                     ! last w-point depth
-            h_oce     = gdepw_0(jpk)
+            IF( rn_bathy > 0.01 ) THEN 
+               IF(lwp) WRITE(numout,*) '         Depth = rn_bathy read in namelist'
+               zdta(:,:) = rn_bathy
+               IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
+                  idta(:,:) = jpkm1
+               ELSE                                                ! z-coordinate (zco or zps): step-like topography
+                  idta(:,:) = jpkm1
+                  DO jk = 1, jpkm1
+                     WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) )   idta(:,:) = jk
+                  END DO
+               ENDIF
+            ELSE
+               IF(lwp) WRITE(numout,*) '         Depth = depthw(jpkm1)'
+               idta(:,:) = jpkm1                            ! before last level
+               zdta(:,:) = gdepw_0(jpk)                     ! last w-point depth
+               h_oce     = gdepw_0(jpk)
+            ENDIF
          ELSE                                         ! bump centered in the basin
             IF(lwp) WRITE(numout,*)
@@ -1106,5 +1120,7 @@
       INTEGER  ::   ios                      ! Local integer output status for namelist read
       REAL(wp) ::   zrmax, ztaper   ! temporary scalars
-      !
+      REAL(wp) ::   zrfact
+      !
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
 
@@ -1115,5 +1131,5 @@
       IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
       !
-      CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
+      CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
       !
       REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
@@ -1169,11 +1185,25 @@
       !                                        ! =============================
       ! use r-value to create hybrid coordinates
+      zenv(:,:) = bathy(:,:)
+      !
+      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
-         END DO
-      END DO
+           IF( bathy(ji,jj) == 0._wp ) THEN
+             iip1 = MIN( ji+1, jpi )
+             ijp1 = MIN( jj+1, jpj )
+             iim1 = MAX( ji-1, 1 )
+             ijm1 = MAX( jj-1, 1 )
+             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
+        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
+               zenv(ji,jj) = rn_sbot_min
+             ENDIF
+           ENDIF
+         END DO
+      END DO
+      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
+      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
       ! 
-      ! Smooth the bathymetry (if required)
+      ! smooth the bathymetry (if required)
       scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
       scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
@@ -1181,78 +1211,68 @@
       jl = 0
       zrmax = 1._wp
-      !                                                     ! ================ !
-      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  !
-         !                                                  ! ================ !
+      !   
+      !     
+      ! set scaling factor used in reducing vertical gradients
+      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
+      !
+      ! initialise temporary evelope depth arrays
+      ztmpi1(:,:) = zenv(:,:)
+      ztmpi2(:,:) = zenv(:,:)
+      ztmpj1(:,:) = zenv(:,:)
+      ztmpj2(:,:) = zenv(:,:)
+      !
+      ! initialise temporary r-value arrays
+      zri(:,:) = 1._wp
+      zrj(:,:) = 1._wp
+      !                                                            ! ================ !
+      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
+         !                                                         ! ================ !
          jl = jl + 1
          zrmax = 0._wp
-         zmsk(:,:) = 0._wp
+         ! we set zrmax from previous r-values (zri and zrj) first
+         ! if set after current r-value calculation (as previously)
+         ! we could exit DO WHILE prematurely before checking r-value
+         ! of current zenv
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
+            END DO
+         END DO
+         zri(:,:) = 0._wp
+         zrj(:,:) = 0._wp
          DO jj = 1, nlcj
             DO ji = 1, nlci
                iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
                ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
-               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
-               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
-               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
-               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
-               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
-               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
-               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
+               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
+                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
+               END IF
+               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
+                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
+               END IF
+               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
+               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
+               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
+               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
             END DO
          END DO
          IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
-         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
-         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp )
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
-                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
-            END DO
-         END DO
          !
-         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
+         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
          !
          DO jj = 1, nlcj
             DO ji = 1, nlci
-               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
-               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
-               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
-               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
-               IF( zmsk(ji,jj) == 1._wp ) THEN
-                  ztmp(ji,jj) =   (                                                                                   &
-             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
-             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
-             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
-             &                    ) / (                                                                               &
-             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
-             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   &
-             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
-             &                        )
-               ENDIF
-            END DO
-         END DO
-         !
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
-               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
-            END DO
-         END DO
-         !
-         ! Apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
-         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp )
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
-               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj)
-            END DO
-         END DO
+               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
+            END DO
+         END DO
+         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
+         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
          !                                                  ! ================ !
       END DO                                                !     End loop     !
       !                                                     ! ================ !
-      !
-      ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions
-      DO ji = nlci+1, jpi 
-         zenv(ji,1:nlcj) = zenv(nlci,1:nlcj)
-      END DO
-      !
-      DO jj = nlcj+1, jpj
-         zenv(:,jj) = zenv(:,nlcj)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
+         END DO
       END DO
       !
@@ -1532,5 +1552,5 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
+      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
