Index: /trunk/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1082)
+++ /trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1083)
@@ -14,4 +14,5 @@
    !!       zgr_zps      : z-coordinate with partial steps
    !!       zgr_sco      : s-coordinate
+   !!       fssig        : sigma coordinate non-dimensional function
    !!---------------------------------------------------------------------
    !! * Modules used
@@ -729,13 +730,466 @@
    END SUBROUTINE zgr_zco
 
-
-
-#  include "domzgr_zps.h90"
-
-
-#if ! defined key_zco
+#if defined key_zco
    !!----------------------------------------------------------------------
-   !!   NOT 'key_zco' :                                    3D gdep. and e3.
+   !!   'key_zco' :                                              "pure" zco (gdep & e3 are 1D arrays)
    !!----------------------------------------------------------------------
+   SUBROUTINE zgr_zps      ! Empty routine
+   END SUBROUTINE zgr_zps
+   SUBROUTINE zgr_sco      ! Empty routine
+   END SUBROUTINE zgr_sco
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                      zco, zps and/or sco available (gedp & e3 are 3D arrays)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE zgr_zps
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zgr_zps  ***
+      !!                     
+      !! ** Purpose :   the depth and vertical scale factor in partial step
+      !!      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(i,j,k)  = fsdep(k)
+      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
+      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
+      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(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(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_0, gdepw_0 and e3._0 are positives
+      !!         - - - - - - -   gdept, gdepw and e3. are positives
+      !!      
+      !!  Reference :
+      !!     Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
+      !!
+      !!  History :
+      !!    8.5  !  02-09 (A. Bozec, G. Madec)  F90: Free form and module
+      !!    9.0  !  02-09 (A. de Miranda)  rigid-lid + islands
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
+      INTEGER  ::   &
+         ik, it            ! temporary integers
+      
+      REAL(wp) ::   &  
+         ze3tp, ze3wp,    &  ! Last ocean level thickness at T- and W-points
+         zdepwp,          &  ! Ajusted ocean depth to avoid too small e3t
+         zdepth,          &  !    "         "
+         zmax, zmin,      &  ! Maximum and minimum depth
+         zdiff               ! temporary scalar
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   &
+         zprt  !    "           "
+
+      LOGICAL ::  ll_print                          ! Allow  control print for debugging
+
+      !!---------------------------------------------------------------------
+      !! OPA8.5, LODYC-IPSL (2002)
+      !!---------------------------------------------------------------------
+      
+      ! Local variable for debugging
+      ll_print=.FALSE.
+!!!   ll_print=.TRUE.
+      
+      ! Initialization of constant
+      zmax = gdepw_0(jpk) + e3t_0(jpk)
+      zmin = gdepw_0(4)
+      
+      ! Ocean depth
+      IF(lwp .AND. ll_print) THEN 
+         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
+
+      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'
+
+
+      ! bathymetry in level (from bathy_meter)
+      ! ===================
+
+      ! initialize mbathy to the maximum ocean level available
+      mbathy(:,:) = jpkm1
+
+      ! storage of land and island's number (zera and negative values) in mbathy
+      DO jj = 1, jpj
+         DO ji= 1, jpi
+            IF( bathy(ji,jj) <= 0. )   mbathy(ji,jj) = INT( bathy(ji,jj) )
+         END DO
+      END DO
+
+      ! bounded value of bathy
+      ! minimum depth == 3 levels
+      ! maximum depth == gdepw_0(jpk)+e3t_0(jpk) 
+      ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk)
+      DO jj = 1, jpj
+         DO ji= 1, jpi
+            IF( bathy(ji,jj) <= 0. ) THEN
+               bathy(ji,jj) = 0.e0
+            ELSE
+               bathy(ji,jj) = MAX( bathy(ji,jj), zmin )
+               bathy(ji,jj) = MIN( bathy(ji,jj), zmax )
+            ENDIF
+         END DO
+      END DO
+
+      ! 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_0 (where
+      ! e3t_0 is the reference level thickness
+      DO jk = jpkm1, 1, -1
+         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1
+            END DO
+         END DO
+      END DO
+
+     ! Scale factors and depth at T- and W-points
+     
+     ! intitialization to the reference z-coordinate
+     DO jk = 1, jpk
+        gdept(:,:,jk) = gdept_0(jk)
+        gdepw(:,:,jk) = gdepw_0(jk)
+        e3t  (:,:,jk) = e3t_0  (jk)
+        e3w  (:,:,jk) = e3w_0  (jk)
+     END DO
+     hdept(:,:) = gdept(:,:,2 )
+     hdepw(:,:) = gdepw(:,:,3 )     
+     
+     ! 
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           ik = mbathy(ji,jj)
+           ! ocean point only
+           IF( ik > 0 ) THEN
+              ! max ocean level case
+              IF( ik == jpkm1 ) THEN
+                 zdepwp = bathy(ji,jj)
+                 ze3tp  = bathy(ji,jj) - gdepw_0(ik)
+                 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
+                 e3t(ji,jj,ik  ) = ze3tp
+                 e3t(ji,jj,ik+1) = ze3tp
+                 e3w(ji,jj,ik  ) = ze3wp
+                 e3w(ji,jj,ik+1) = ze3tp
+                 gdepw(ji,jj,ik+1) = zdepwp
+                 gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
+                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
+                 ! standard case
+              ELSE
+!!alex
+                 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN
+                    gdepw(ji,jj,ik+1) = bathy(ji,jj)
+                 ELSE
+!!alex ctl          write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1)
+                    gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
+                 ENDIF
+!!Alex
+!!Alex           gdepw(ji,jj,ik+1) = bathy(ji,jj)
+!!bug gm  verifier les gdepw_0
+                 gdept(ji,jj,ik  ) =  gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) )   &
+                    &                             * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
+                    &                             / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
+                 e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
+                    &                      / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
+                 e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   &
+                    &                * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
+                 !       ... on ik+1
+                 e3w(ji,jj,ik+1) = e3t(ji,jj,ik)
+                 e3t(ji,jj,ik+1) = e3t(ji,jj,ik)
+                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
+              ENDIF
+           ENDIF
+        END DO
+     END DO
+
+     it = 0
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           ik = mbathy(ji,jj)
+           ! ocean point only
+           IF( ik > 0 ) THEN
+              ! bathymetry output
+              hdept(ji,jj) = gdept(ji,jj,ik  )
+              hdepw(ji,jj) = gdepw(ji,jj,ik+1)
+              e3tp (ji,jj) = e3t(ji,jj,ik  )
+              e3wp (ji,jj) = e3w(ji,jj,ik  )
+              ! test
+              zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
+              IF( zdiff <= 0. .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= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1),   &
+                                 ' zdiff   = ', zdiff
+                 WRITE(numout,*) ' e3tp    = ', e3t(ji,jj,ik  ), ' e3wp    = ', e3w(ji,jj,ik  )
+              ENDIF
+           ENDIF
+        END DO
+      END DO
+
+      ! Scale factors and depth at U-, V-, UW and VW-points
+
+      ! initialisation to z-scale factors
+      DO jk = 1, jpk
+         e3u (:,:,jk)  = e3t_0(jk)
+         e3v (:,:,jk)  = e3t_0(jk)
+         e3uw(:,:,jk)  = e3w_0(jk)
+         e3vw(:,:,jk)  = e3w_0(jk)
+      END DO
+
+     ! Computed as the minimum of neighbooring scale factors
+     DO jk = 1,jpk
+        DO jj = 1, jpjm1
+           DO ji = 1, fs_jpim1   ! vector opt.
+              e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
+              e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
+              e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
+              e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
+           END DO
+        END DO
+     END DO
+     
+     ! Boundary conditions
+     CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. )
+     CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. )
+     
+     ! set to z-scale factor if zero (i.e. along closed boundaries)
+     DO jk = 1, jpk
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk)
+              IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk)
+              IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk)
+              IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk)
+           END DO
+        END DO
+     END DO
+     
+     ! Scale factor at F-point
+     
+     ! initialisation to z-scale factors
+     DO jk = 1, jpk
+        e3f (:,:,jk) = e3t_0(jk)
+     END DO
+     
+     ! Computed as the minimum of neighbooring V-scale factors
+     DO jk = 1, jpk
+        DO jj = 1, jpjm1
+           DO ji = 1, fs_jpim1   ! vector opt.
+              e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
+           END DO
+        END DO
+     END DO
+     ! Boundary conditions
+     CALL lbc_lnk( e3f, 'F', 1. )
+     
+     ! set to z-scale factor if zero (i.e. along closed boundaries)
+     DO jk = 1, jpk
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk)
+           END DO
+        END DO
+     END DO
+!!bug ? gm:  must be a do loop with mj0,mj1
+
+     ! we duplicate factor scales for jj = 1 and jj = 2
+     e3t(:,mj0(1),:) = e3t(:,mj0(2),:) 
+     e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
+     e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
+     e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
+     e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
+
+
+     
+     ! Compute gdep3w (vertical sum of e3w)
+     
+     gdep3w   (:,:,1) = 0.5 * e3w (:,:,1)
+     DO jk = 2, jpk
+        gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w (:,:,jk) 
+     END DO
+          
+     ! Control print
+ 9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ')
+ 9610 FORMAT(10x,i4,4f9.2)
+      IF(lwp .AND. ll_print) THEN
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik = MAX(mbathy(ji,jj),1)
+               zprt(ji,jj) = e3t(ji,jj,ik)
+            END DO
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3t(mbathy)'
+         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik = MAX(mbathy(ji,jj),1)
+               zprt(ji,jj) = e3w(ji,jj,ik)
+            END DO
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3w(mbathy)'
+         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik = MAX(mbathy(ji,jj),1)
+               zprt(ji,jj) = e3u(ji,jj,ik)
+            END DO
+         END DO
+
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3u(mbathy)'
+         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik = MAX(mbathy(ji,jj),1)
+               zprt(ji,jj) = e3v(ji,jj,ik)
+            END DO
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3v(mbathy)'
+         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik = MAX(mbathy(ji,jj),1) 
+               zprt(ji,jj) = e3f(ji,jj,ik)
+            END DO
+         END DO
+
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3f(mbathy)'
+         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik =  MAX(mbathy(ji,jj),1)
+               zprt(ji,jj) = gdep3w(ji,jj,ik)
+            END DO
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr gdep3w(mbathy)'
+         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+
+      ENDIF  
+
+
+      DO jk = 1,jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi 
+               IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN
+                  IF(lwp) THEN
+                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
+                     WRITE(numout,*) ' =========           --------------- '
+                     WRITE(numout,*)
+                  ENDIF
+                  STOP 'domzgr.psteps'
+               ENDIF
+               IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN
+                  IF(lwp) THEN
+                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
+                     WRITE(numout,*) ' =========        ------------------ '
+                     WRITE(numout,*)
+                  ENDIF
+                  STOP 'domzgr.psteps'
+               ENDIF
+            END DO
+         END DO
+      END DO  
+
+   IF(lwp) THEN
+      WRITE(numout,*) ' e3t lev 21 '
+      CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3w lev 21  '
+      CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3u lev 21  '
+      CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3v lev 21  '
+      CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3f lev 21  '
+      CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3t lev 22 '
+      CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3w lev 22  '
+      CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3u lev 22  '
+      CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3v lev 22  '
+      CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
+      WRITE(numout,*) ' e3f lev 22  '
+      CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
+   ENDIF
+
+      ! ===========
+      ! Zoom domain 
+      ! ===========
+
+      IF( lzoom )   CALL zgr_bat_zoom
+
+      ! ================
+      ! Bathymetry check
+      ! ================
+
+      CALL zgr_bat_ctl
+
+   END SUBROUTINE zgr_zps
+
+
+   FUNCTION fssig( pk ) RESULT( pf )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_init  ***
+      !!       
+      !! ** Purpose :   provide the analytical function in s-coordinate
+      !!          
+      !! ** Method  :   the function provide the non-dimensional position of
+      !!                T and W (i.e. between 0 and 1)
+      !!                T-points at integer values (between 1 and jpk)
+      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+      !!
+      !! Reference  :   ???
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   pk   ! continuous "k" coordinate
+      REAL(wp)                ::   pf   ! sigma value
+      !!----------------------------------------------------------------------
+      !
+      pf =   (   TANH( theta * ( -(pk-0.5) / REAL(jpkm1) + thetb )  )      &
+         &     - TANH( thetb * theta                                )  )   &
+         & * (   COSH( theta                           )                   &
+         &     + COSH( theta * ( 2.e0 * thetb - 1.e0 ) )  )                &
+         & / ( 2.e0 * SINH( theta ) )
+      !
+   END FUNCTION fssig
+
 
    SUBROUTINE zgr_sco
@@ -758,9 +1212,9 @@
       !!            hbatf = mi( mj( hbatt ) )
       !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
-      !!         function and its derivative given as statement function.
-      !!            gsigt(k) = fssig (k+0.5)
-      !!            gsigw(k) = fssig (k    )
-      !!            esigt(k) = fsdsig(k+0.5)
-      !!            esigw(k) = fsdsig(k    )
+      !!         function and its derivative given as function.
+      !!            gsigt(k) = fssig (k    )
+      !!            gsigw(k) = fssig (k-0.5)
+      !!            esigt(k) = fsdsig(k    )
+      !!            esigw(k) = fsdsig(k-0.5)
       !!      This routine is given as an example, it must be modified
       !!      following the user s desiderata. nevertheless, the output as
@@ -781,10 +1235,6 @@
       !! * Local declarations
       INTEGER  ::   ji, jj, jk, jl
-      REAL(wp) ::   fssig, fsdsig, pfkk
-
       INTEGER  ::   iip1, ijp1, iim1, ijm1
-      REAL(wp) ::   &
-         fssigt, fssigw, zcoeft, zcoefw,   &
-         zrmax, ztaper
+      REAL(wp) ::   zrmax, ztaper
 
       REAL(wp), DIMENSION(jpi,jpj) ::   &
@@ -792,15 +1242,4 @@
 
       NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
-      !!----------------------------------------------------------------------
-
-      ! a. Sigma stretching coefficient
-       fssigt(pfkk) = (tanh(theta*(-(pfkk-0.5)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &
-                    * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))
-       fssigw(pfkk) = (tanh(theta*(-(pfkk-1.0)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &
-                    * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))
-
-      ! b. Vertical derivative of sigma stretching coefficient
- 
-!!bug fsdsig(pfkk)= put here the analytical derivative of fssigt and w ...
       !!----------------------------------------------------------------------
 
@@ -1043,16 +1482,10 @@
       ! 2.1 Coefficients for model level depth at w- and t-levels
 
-!!bug gm : change the def of statment function....
       DO jk = 1, jpk
-        gsigw(jk) = -fssigw(FLOAT(jk))
-        gsigt(jk) = -fssigt(FLOAT(jk))
+        gsigw(jk) = -fssig( FLOAT(jk)-0.5 )
+        gsigt(jk) = -fssig( FLOAT(jk)     )
       END DO
 
       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
-
-!!org DO jk = 1, jpk
-!!org    gsigw(jk) = -fssig ( FLOAT(jk)    )
-!!org    gsigt(jk) = -fssig ( FLOAT(jk)+0.5)
-!!org END DO
 
       ! 2.2 Coefficients for vertical scale factors at w-, t- levels
@@ -1068,9 +1501,4 @@
         esigw(1)=esigw(2)
         esigt(jpk)=esigt(jpkm1)
-
-!!org DO jk = 1, jpk
-!!org    esigw(jk)=fsdsig( FLOAT(jk))
-!!org    esigt(jk)=fsdsig( FLOAT(jk)+0.5)
-!!org END DO
 
       ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors
@@ -1205,10 +1633,4 @@
    END SUBROUTINE zgr_sco
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                      Empty routine
-   !!----------------------------------------------------------------------
-   SUBROUTINE zgr_sco
-   END SUBROUTINE zgr_sco
 #endif
 
Index: unk/NEMO/OPA_SRC/DOM/domzgr_zps.h90
===================================================================
--- /trunk/NEMO/OPA_SRC/DOM/domzgr_zps.h90	(revision 1082)
+++ 	(revision )
@@ -1,443 +1,0 @@
-   !!----------------------------------------------------------------------
-   !!                      ***  domzgr_zps.h90  ***
-   !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-
-#if ! defined key_zco
-   !!----------------------------------------------------------------------
-   !!   NOT 'key_zco' :                                    3D gdep. and e3.
-   !!----------------------------------------------------------------------
-
-   SUBROUTINE zgr_zps
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zgr_zps  ***
-      !!                     
-      !! ** Purpose :   the depth and vertical scale factor in partial step
-      !!      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(i,j,k)  = fsdep(k)
-      !!                       e3w(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
-      !!              t-level: gdept(i,j,k)  = fsdep(k+0.5)
-      !!                       e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(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(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_0, gdepw_0 and e3._0 are positives
-      !!         - - - - - - -   gdept, gdepw and e3. are positives
-      !!      
-      !!  Reference :
-      !!     Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
-      !!
-      !!  History :
-      !!    8.5  !  02-09 (A. Bozec, G. Madec)  F90: Free form and module
-      !!    9.0  !  02-09 (A. de Miranda)  rigid-lid + islands
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      INTEGER  ::   &
-         ik, it            ! temporary integers
-      
-      REAL(wp) ::   &  
-         ze3tp, ze3wp,    &  ! Last ocean level thickness at T- and W-points
-         zdepwp,          &  ! Ajusted ocean depth to avoid too small e3t
-         zdepth,          &  !    "         "
-         zmax, zmin,      &  ! Maximum and minimum depth
-         zdiff               ! temporary scalar
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zprt  !    "           "
-
-      LOGICAL ::  ll_print                          ! Allow  control print for debugging
-
-      !!---------------------------------------------------------------------
-      !! OPA8.5, LODYC-IPSL (2002)
-      !!---------------------------------------------------------------------
-      
-      ! Local variable for debugging
-      ll_print=.FALSE.
-!!!   ll_print=.TRUE.
-      
-      ! Initialization of constant
-      zmax = gdepw_0(jpk) + e3t_0(jpk)
-      zmin = gdepw_0(4)
-      
-      ! Ocean depth
-      IF(lwp .AND. ll_print) THEN 
-         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
-
-      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'
-
-
-      ! bathymetry in level (from bathy_meter)
-      ! ===================
-
-      ! initialize mbathy to the maximum ocean level available
-      mbathy(:,:) = jpkm1
-
-      ! storage of land and island's number (zera and negative values) in mbathy
-      DO jj = 1, jpj
-         DO ji= 1, jpi
-            IF( bathy(ji,jj) <= 0. )   mbathy(ji,jj) = INT( bathy(ji,jj) )
-         END DO
-      END DO
-
-      ! bounded value of bathy
-      ! minimum depth == 3 levels
-      ! maximum depth == gdepw_0(jpk)+e3t_0(jpk) 
-      ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk)
-      DO jj = 1, jpj
-         DO ji= 1, jpi
-            IF( bathy(ji,jj) <= 0. ) THEN
-               bathy(ji,jj) = 0.e0
-            ELSE
-               bathy(ji,jj) = MAX( bathy(ji,jj), zmin )
-               bathy(ji,jj) = MIN( bathy(ji,jj), zmax )
-            ENDIF
-         END DO
-      END DO
-
-      ! 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_0 (where
-      ! e3t_0 is the reference level thickness
-      DO jk = jpkm1, 1, -1
-         zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth )   mbathy(ji,jj) = jk-1
-            END DO
-         END DO
-      END DO
-
-     ! Scale factors and depth at T- and W-points
-     
-     ! intitialization to the reference z-coordinate
-     DO jk = 1, jpk
-        gdept(:,:,jk) = gdept_0(jk)
-        gdepw(:,:,jk) = gdepw_0(jk)
-        e3t  (:,:,jk) = e3t_0  (jk)
-        e3w  (:,:,jk) = e3w_0  (jk)
-     END DO
-     hdept(:,:) = gdept(:,:,2 )
-     hdepw(:,:) = gdepw(:,:,3 )     
-     
-     ! 
-     DO jj = 1, jpj
-        DO ji = 1, jpi
-           ik = mbathy(ji,jj)
-           ! ocean point only
-           IF( ik > 0 ) THEN
-              ! max ocean level case
-              IF( ik == jpkm1 ) THEN
-                 zdepwp = bathy(ji,jj)
-                 ze3tp  = bathy(ji,jj) - gdepw_0(ik)
-                 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
-                 e3t(ji,jj,ik  ) = ze3tp
-                 e3t(ji,jj,ik+1) = ze3tp
-                 e3w(ji,jj,ik  ) = ze3wp
-                 e3w(ji,jj,ik+1) = ze3tp
-                 gdepw(ji,jj,ik+1) = zdepwp
-                 gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
-                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
-                 ! standard case
-              ELSE
-!!alex
-                 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN
-                    gdepw(ji,jj,ik+1) = bathy(ji,jj)
-                 ELSE
-!!alex ctl          write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1)
-                    gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
-                 ENDIF
-!!Alex
-!!Alex           gdepw(ji,jj,ik+1) = bathy(ji,jj)
-!!bug gm  verifier les gdepw_0
-                 gdept(ji,jj,ik  ) =  gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) )   &
-                    &                             * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
-                    &                             / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
-                 e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
-                    &                      / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
-                 e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   &
-                    &                * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
-                 !       ... on ik+1
-                 e3w(ji,jj,ik+1) = e3t(ji,jj,ik)
-                 e3t(ji,jj,ik+1) = e3t(ji,jj,ik)
-                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
-              ENDIF
-           ENDIF
-        END DO
-     END DO
-
-     it = 0
-     DO jj = 1, jpj
-        DO ji = 1, jpi
-           ik = mbathy(ji,jj)
-           ! ocean point only
-           IF( ik > 0 ) THEN
-              ! bathymetry output
-              hdept(ji,jj) = gdept(ji,jj,ik  )
-              hdepw(ji,jj) = gdepw(ji,jj,ik+1)
-              e3tp (ji,jj) = e3t(ji,jj,ik  )
-              e3wp (ji,jj) = e3w(ji,jj,ik  )
-              ! test
-              zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
-              IF( zdiff <= 0. .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= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1),   &
-                                 ' zdiff   = ', zdiff
-                 WRITE(numout,*) ' e3tp    = ', e3t(ji,jj,ik  ), ' e3wp    = ', e3w(ji,jj,ik  )
-              ENDIF
-           ENDIF
-        END DO
-      END DO
-
-      ! Scale factors and depth at U-, V-, UW and VW-points
-
-      ! initialisation to z-scale factors
-      DO jk = 1, jpk
-         e3u (:,:,jk)  = e3t_0(jk)
-         e3v (:,:,jk)  = e3t_0(jk)
-         e3uw(:,:,jk)  = e3w_0(jk)
-         e3vw(:,:,jk)  = e3w_0(jk)
-      END DO
-
-     ! Computed as the minimum of neighbooring scale factors
-     DO jk = 1,jpk
-        DO jj = 1, jpjm1
-           DO ji = 1, fs_jpim1   ! vector opt.
-              e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
-              e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
-              e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
-              e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
-           END DO
-        END DO
-     END DO
-     
-     ! Boundary conditions
-     CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. )
-     CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. )
-     
-     ! set to z-scale factor if zero (i.e. along closed boundaries)
-     DO jk = 1, jpk
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk)
-              IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk)
-              IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk)
-              IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk)
-           END DO
-        END DO
-     END DO
-     
-     ! Scale factor at F-point
-     
-     ! initialisation to z-scale factors
-     DO jk = 1, jpk
-        e3f (:,:,jk) = e3t_0(jk)
-     END DO
-     
-     ! Computed as the minimum of neighbooring V-scale factors
-     DO jk = 1, jpk
-        DO jj = 1, jpjm1
-           DO ji = 1, fs_jpim1   ! vector opt.
-              e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
-           END DO
-        END DO
-     END DO
-     ! Boundary conditions
-     CALL lbc_lnk( e3f, 'F', 1. )
-     
-     ! set to z-scale factor if zero (i.e. along closed boundaries)
-     DO jk = 1, jpk
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk)
-           END DO
-        END DO
-     END DO
-!!bug ? gm:  must be a do loop with mj0,mj1
-
-     ! we duplicate factor scales for jj = 1 and jj = 2
-     e3t(:,mj0(1),:) = e3t(:,mj0(2),:) 
-     e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
-     e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
-     e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
-     e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
-
-
-     
-     ! Compute gdep3w (vertical sum of e3w)
-     
-     gdep3w   (:,:,1) = 0.5 * e3w (:,:,1)
-     DO jk = 2, jpk
-        gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w (:,:,jk) 
-     END DO
-          
-     ! Control print
- 9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ')
- 9610 FORMAT(10x,i4,4f9.2)
-      IF(lwp .AND. ll_print) THEN
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3t(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3t(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3w(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3w(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3u(ji,jj,ik)
-            END DO
-         END DO
-
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3u(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3v(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3v(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1) 
-               zprt(ji,jj) = e3f(ji,jj,ik)
-            END DO
-         END DO
-
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3f(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik =  MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = gdep3w(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr gdep3w(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-
-      ENDIF  
-
-
-      DO jk = 1,jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi 
-               IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN
-                  IF(lwp) THEN
-                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
-                     WRITE(numout,*) ' =========           --------------- '
-                     WRITE(numout,*)
-                  ENDIF
-                  STOP 'domzgr.psteps'
-               ENDIF
-               IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN
-                  IF(lwp) THEN
-                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
-                     WRITE(numout,*) ' =========        ------------------ '
-                     WRITE(numout,*)
-                  ENDIF
-                  STOP 'domzgr.psteps'
-               ENDIF
-            END DO
-         END DO
-      END DO  
-
-   IF(lwp) THEN
-      WRITE(numout,*) ' e3t lev 21 '
-      CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3w lev 21  '
-      CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3u lev 21  '
-      CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3v lev 21  '
-      CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3f lev 21  '
-      CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3t lev 22 '
-      CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3w lev 22  '
-      CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3u lev 22  '
-      CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3v lev 22  '
-      CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3f lev 22  '
-      CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-   ENDIF
-
-      ! ===========
-      ! Zoom domain 
-      ! ===========
-
-      IF( lzoom )   CALL zgr_bat_zoom
-
-      ! ================
-      ! Bathymetry check
-      ! ================
-
-      CALL zgr_bat_ctl
-
-   END SUBROUTINE zgr_zps
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                      Empty routine
-   !!----------------------------------------------------------------------
-   SUBROUTINE zgr_zps
-   END SUBROUTINE zgr_zps
-#endif
