Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3226)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3432)
@@ -54,5 +54,5 @@
    !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
    REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates
-
+   PUBLIC rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb,rn_hc
    !! * Control permutation of array indices
 #  include "oce_ftrans.h90"
@@ -82,5 +82,5 @@
       !!                   ln_zco=T   z-coordinate   
       !!                   ln_zps=T   z-coordinate with partial steps
-      !!                   ln_sco=T   s-coordinate 
+      !!                   ln_zco=T   s-coordinate 
       !!
       !! ** Action  :   define gdep., e3., mbathy and bathy
@@ -395,6 +395,11 @@
          mbathy(:,:) = 0                                   ! set to zero extra halo points
          bathy (:,:) = 0._wp                               ! (require for mpp case)
+#if defined key_mpp_rkpart
+         DO jj = nldj, nlcj                                   ! interior values
+            DO ji = nldi, nlci
+#else
          DO jj = 1, nlcj                                   ! interior values
             DO ji = 1, nlci
+#endif
                mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
                bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
@@ -593,6 +598,7 @@
       !!              - update bathy : meter bathymetry (in meters)
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zbathy => wrk_2d_1
+      USE mapcomm_mod, ONLY:   trimmed, nidx,eidx,sidx,widx
+      USE wrk_nemo,    ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo,    ONLY:   zbathy => wrk_2d_1
       !!
       INTEGER ::   ji, jj, jl                    ! dummy loop indices
@@ -636,4 +642,5 @@
          IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
       ENDIF
+
       IF( lk_mpp ) THEN
          zbathy(:,:) = FLOAT( mbathy(:,:) )
@@ -646,8 +653,8 @@
          IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
          IF( lk_mpp ) THEN
-            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
+            IF( (nbondi == -1 .OR. nbondi == 2) .AND. (.NOT. trimmed(widx,narea) ) ) THEN
                IF( jperio /= 1 )   mbathy(1,:) = 0
             ENDIF
-            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
+            IF( (nbondi == 1 .OR. nbondi == 2) .AND. (.NOT. trimmed(eidx,narea) ) ) THEN
                IF( jperio /= 1 )   mbathy(nlci,:) = 0
             ENDIF
@@ -1188,5 +1195,5 @@
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
       USE wrk_nemo, ONLY:   zenv => wrk_2d_1 , ztmp => wrk_2d_2 , zmsk  => wrk_2d_3
-      USE wrk_nemo, ONLY:   zri  => wrk_2d_4 , zrj  => wrk_2d_5 , zhbat => wrk_2d_6
+      USE wrk_nemo, ONLY:   ztmp2 => wrk_2d_4 , zhbat => wrk_2d_5
       USE wrk_nemo, ONLY:   gsigw3  => wrk_3d_1
       USE wrk_nemo, ONLY:   gsigt3  => wrk_3d_2
@@ -1199,4 +1206,7 @@
       USE wrk_nemo, ONLY:   esigwu3 => wrk_3d_9
       USE wrk_nemo, ONLY:   esigwv3 => wrk_3d_10
+      USE mapcomm_mod, ONLY: trimmed, cyclic_bc
+      USE mapcomm_mod, ONLY: nidx, eidx, sidx, widx
+!      USE arpdebugging, ONLY: dump_array
       !! DCSE_NEMO: wrk_nemo module variables renamed, need additional directives
 !FTRANS gsigw3 :I :I :z
@@ -1213,5 +1223,6 @@
       INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
       INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
-      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars
+      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper, zri, zrj   ! temporary scalars
+      REAL(wp), PARAMETER :: TOL_ZERO = 1.0E-20_wp ! Any value less than this assumed zero
       !
 
@@ -1219,5 +1230,5 @@
       !!----------------------------------------------------------------------
 
-      IF( wrk_in_use(2, 1,2,3,4,5,6) .OR. wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10) ) THEN
+      IF( wrk_in_use(2, 1,2,3,4,5) .OR. wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10) ) THEN
          CALL ctl_stop('zgr_sco: ERROR - requested workspace arrays unavailable')   ;   RETURN
       ENDIF
@@ -1269,4 +1280,6 @@
          END DO
       END DO
+
+      CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. )
       ! 
       ! Smooth the bathymetry (if required)
@@ -1282,30 +1295,39 @@
          zrmax = 0._wp
          zmsk(:,:) = 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
-            END DO
-         END DO
+               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last row  (jj=nclj+1 to jpj)
+               zri = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
+               zrj = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
+               zrmax = MAX( zrmax, zri, zrj )
+               IF( zri > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
+               IF( zri > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
+               IF( zrj > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
+               IF( zrj > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
+            END DO
+         END DO
+
+         ! lateral boundary condition on zmsk: retain any 1's along closed 
+         ! boundary (use of lzero flag to lbc_lnk)
+         CALL lbc_lnk( zmsk, 'T', 1._wp, lzero=.FALSE. )
+
          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=',I5,' rmax=',F8.4,' nb of pt= ',I8)") &
+                                                         jl, zrmax, INT( SUM(zmsk(:,:) ) )
          !
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
+!!$         IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN
+!!$            CALL dump_array(jl, 'zenv_before', zenv, withHalos=.TRUE.)
+!!$            CALL dump_array(jl, 'ztmp_before', ztmp, withHalos=.TRUE.)
+!!$            CALL dump_array(jl, 'zmsk_before', zmsk, withHalos=.TRUE.)
+!!$         END IF
+
+         ! Copy current surface before next smoothing iteration 
+         ztmp(:,:) = zenv(:,:)
+
+         DO jj = nldj, nlcj
+            DO ji = nldi, nlci
                iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
                ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
@@ -1326,17 +1348,24 @@
          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) )
+
+         ! Need to update halos of ztmp here but do not zero halos on closed 
+         ! boundaries
+         CALL lbc_lnk( ztmp, 'T', 1._wp, lzero=.FALSE.)
+
+         DO jj = 1,nlcj
+            DO ji = 1,nlci
+               IF( zmsk(ji,jj) >= 1._wp-TOL_ZERO ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
             END DO
          END DO
          !
-         ! Apply lateral boundary condition   CAUTION: kept 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
+         ! Apply lateral boundary condition but do not zero on closed boundaries
+         CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. )
+
+!!$         IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN
+!!$            CALL dump_array(jl, 'zenv', zenv, withHalos=.TRUE.)
+!!$            CALL dump_array(jl, 'ztmp', ztmp, withHalos=.TRUE.)
+!!$            CALL dump_array(jl, 'zmsk', zmsk, withHalos=.TRUE.)
+!!$         END IF
+
          !                                                  ! ================ !
       END DO                                                !     End loop     !
@@ -1365,4 +1394,6 @@
          ENDIF
       ENDIF
+
+!      CALL dump_array(0, 'hbatt', hbatt, withHalos=.FALSE.)
 
       !                                        ! ==============================
@@ -1576,4 +1607,12 @@
       !
 !!    H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
+      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0
+      where (e3u   (:,:,:).eq.0.0)  e3u(:,:,:) = 1.0
+      where (e3v   (:,:,:).eq.0.0)  e3v(:,:,:) = 1.0
+      where (e3f   (:,:,:).eq.0.0)  e3f(:,:,:) = 1.0
+      where (e3w   (:,:,:).eq.0.0)  e3w(:,:,:) = 1.0
+      where (e3uw  (:,:,:).eq.0.0)  e3uw(:,:,:) = 1.0
+      where (e3vw  (:,:,:).eq.0.0)  e3vw(:,:,:) = 1.0
+
 
       fsdept(:,:,:) = gdept (:,:,:)
@@ -1597,4 +1636,5 @@
          END DO
       END DO
+
       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
          &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
@@ -1642,14 +1682,15 @@
             END DO
          END DO
-         DO jj = mj0(74), mj1(74)
-            DO ji = mi0(100), mi1(100)
-               WRITE(numout,*)
-               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
-               WRITE(numout,*) ' ~~~~~~  --------------------'
-               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
-               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
-                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
-            END DO
-         END DO
+!!$ ARPDBG - out of bounds if jpj < 74 or jpi < 100, e.g. default GYRE
+!!$         DO jj = mj0(74), mj1(74)
+!!$            DO ji = mi0(100), mi1(100)
+!!$               WRITE(numout,*)
+!!$               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
+!!$               WRITE(numout,*) ' ~~~~~~  --------------------'
+!!$               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
+!!$               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
+!!$                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
+!!$            END DO
+!!$         END DO
       ENDIF
 
@@ -1668,5 +1709,5 @@
                   CALL ctl_stop( ctmp1 )
                ENDIF
-               IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
+               IF( gdepw_1(ji,jj,jk) < 0._wp .OR. gdept_1(ji,jj,jk) < 0._wp ) THEN
                   WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
                   CALL ctl_stop( ctmp1 )
@@ -1677,5 +1718,5 @@
 !!gm bug    #endif
       !
-      IF( wrk_not_released(2, 1,2,3,4,5,6) .OR. wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10) )  &
+      IF( wrk_not_released(2, 1,2,3,4,5) .OR. wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10) )  &
         &  CALL ctl_stop('dom:zgr_sco: failed to release workspace arrays')
       !
