Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3925)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3926)
@@ -169,5 +169,5 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zgr_z  ***
-      !!                   
+      !!                    
       !! ** Purpose :   set the depth of model levels and the resulting 
       !!      vertical scale factors.
@@ -639,4 +639,5 @@
          END DO
       END DO
+      IF( lk_mpp )   CALL mpp_sum( icompt )
       IF( icompt == 0 ) THEN
          IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
@@ -1101,6 +1102,9 @@
       INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
       REAL(wp) ::   zrmax, ztaper   ! temporary scalars
-      !
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
+      REAL(wp) ::   zrfact   ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
+
+      !
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, zri, zrj, zhbat
 
       NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
@@ -1110,6 +1114,7 @@
       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,      ztmpi1, ztmpi2, ztmpj1, ztmpj2         )
+      CALL wrk_alloc( jpi, jpj,      zenv, zri, zrj, zhbat     )
+     !
       REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters
       READ  ( numnam, namzgr_sco )
@@ -1158,9 +1163,11 @@
       !                                        ! =============================
       ! use r-value to create hybrid coordinates
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
-         END DO
-      END DO
+!     DO jj = 1, jpj
+!        DO ji = 1, jpi
+!           zenv(ji,jj) = MAX( bathy(ji,jj), 0._wp )
+!        END DO
+!     END DO
+!     CALL lbc_lnk( zenv, 'T', 1._wp )
+      zenv(:,:) = bathy(:,:)
       ! 
       ! Smooth the bathymetry (if required)
@@ -1170,82 +1177,73 @@
       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 abd 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
+               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
             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
+         CALL lbc_lnk( zenv, 'T', 1._wp )
          !                                                  ! ================ !
       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)
-      END DO
+!     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 values
+!        END DO
+!     END DO
       !
       ! Envelope bathymetry saved in hbatt
       hbatt(:,:) = zenv(:,:) 
+
       IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
          CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
@@ -1493,6 +1491,5 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
-      !
+      CALL wrk_dealloc( jpi, jpj,      zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat                           )      !
       IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
       !
