Index: /trunk/CONFIG/GYRE/EXP00/namelist
===================================================================
--- /trunk/CONFIG/GYRE/EXP00/namelist	(revision 1347)
+++ /trunk/CONFIG/GYRE/EXP00/namelist	(revision 1348)
@@ -64,4 +64,7 @@
    thetb       =    0.75   !  bottom control parameter  (0<=thetb<= 1)
    r_max       =    0.15   !  maximum cut-off r-value allowed (0<r_max<1)
+   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
+   bb          =    0.8    !  stretching with s-sigma
+   hc          =  150.0    !  critical depth with s-sigma
 /
 !-----------------------------------------------------------------------
Index: /trunk/CONFIG/GYRE_LOBSTER/EXP00/namelist
===================================================================
--- /trunk/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 1347)
+++ /trunk/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 1348)
@@ -64,4 +64,7 @@
    thetb       =    0.75   !  bottom control parameter  (0<=thetb<= 1)
    r_max       =    0.15   !  maximum cut-off r-value allowed (0<r_max<1)
+   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
+   bb          =    0.8    !  stretching with s-sigma
+   hc          =  150.0    !  critical depth with s-sigma
 /
 !-----------------------------------------------------------------------
Index: /trunk/CONFIG/ORCA2_LIM/EXP00/1_namelist
===================================================================
--- /trunk/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 1347)
+++ /trunk/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 1348)
@@ -64,4 +64,7 @@
    thetb       =    0.75   !  bottom control parameter  (0<=thetb<= 1)
    r_max       =    0.15   !  maximum cut-off r-value allowed (0<r_max<1)
+   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
+   bb          =    0.8    !  stretching with s-sigma
+   hc          =  150.0    !  critical depth with s-sigma
 /
 !-----------------------------------------------------------------------
Index: /trunk/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /trunk/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 1347)
+++ /trunk/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 1348)
@@ -64,4 +64,7 @@
    thetb       =    0.75   !  bottom control parameter  (0<=thetb<= 1)
    r_max       =    0.15   !  maximum cut-off r-value allowed (0<r_max<1)
+   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
+   bb          =    0.8    !  stretching with s-sigma
+   hc          =  150.0    !  critical depth with s-sigma 
 /
 !-----------------------------------------------------------------------
Index: /trunk/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist
===================================================================
--- /trunk/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 1347)
+++ /trunk/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 1348)
@@ -64,4 +64,7 @@
    thetb       =    0.75   !  bottom control parameter  (0<=thetb<= 1)
    r_max       =    0.15   !  maximum cut-off r-value allowed (0<r_max<1)
+   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
+   bb          =    0.8    !  stretching with s-sigma
+   hc          =  150.0    !  critical depth with s-sigma
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1347)
+++ /trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1348)
@@ -112,5 +112,24 @@
       IF( ln_zps )   CALL zgr_zps      ! Partial step z-coordinate
       IF( ln_sco )   CALL zgr_sco      ! s-coordinate or hybrid z-s coordinate
-      !
+
+!!bug gm control print:
+      IF( nprint == 1 .AND. lwp )   THEN
+         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
+         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
+            &                   ' w ',   MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ),  &
+            &                   ' u ',   MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ),  &
+            &                   ' uw',   MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)),   &
+            &                   ' w ',   MINVAL( fse3w(:,:,:) )
+
+         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
+            &                   ' w ',   MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ),  &
+            &                   ' u ',   MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ),  &
+            &                   ' uw',   MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)),   &
+            &                   ' w ',   MAXVAL( fse3w(:,:,:) )
+      ENDIF
+!!!bug gm
+
    END SUBROUTINE dom_zgr
 
@@ -264,6 +283,6 @@
       INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
       INTEGER  ::   inum                      ! temporary logical unit
+      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
       INTEGER  ::   ii0, ii1, ij0, ij1        ! indices
-      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
       REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics 
       REAL(wp) ::   zi     , zj     , zh      ! temporary scalars
@@ -394,5 +413,5 @@
             CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
             CALL iom_close (inum)
-           !                                                ! =====================
+            !                                                ! =====================
            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
               !                                             ! =====================
@@ -422,6 +441,6 @@
               !
            ENDIF
-
-         ENDIF
+            !
+        ENDIF
          !                                            ! =============== !
       ELSE                                            !      error      !
@@ -996,4 +1015,34 @@
 
 
+   FUNCTION fssig1( pk1, bb ) RESULT( pf1 )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_init  ***
+      !!
+      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
+      !!
+      !! ** Method  :   the function provides 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   ) ::   pk1   ! continuous "k" coordinate
+      REAL(wp), INTENT(in   ) ::   bb    ! Stretching coefficient
+      REAL(wp)                ::   pf1   ! sigma value
+
+      !!----------------------------------------------------------------------
+      !
+      IF ( theta == 0 ) then   !uniform sigma
+         pf1 = -(pk1-0.5)/REAL(jpkm1)
+      ELSE    ! stretched sigma
+         pf1 =   (1.0-bb) * (sinh( theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(theta) + &
+                 bb * ( (tanh( theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*theta) ) / &
+                 (2*tanh(0.5*theta) ) )
+      ENDIF
+
+   END FUNCTION fssig1
+
+
    SUBROUTINE zgr_sco
       !!----------------------------------------------------------------------
@@ -1033,6 +1082,12 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zenv, ztmp, zmsk    ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zri , zrj , zhbat   !  -     -
-      !!
-      NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
+
+      LOGICAL :: ln_s_sigma = .false. !use hybrid s_sigma coordinates & stretching function fssig1,used with ln_sco = .true.
+      REAL(wp) :: bb = 0.8   ! stretching parameter for song and haidvogel stretching, bb=0; top only, bb =1; top and bottom
+      REAL(wp) :: hc = 150   ! Critical depth for s-sigma coordinates
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   gsigw3, gsigt3, gsi3w3, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3
+      !!
+      NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max, ln_s_sigma, bb, hc
       !!----------------------------------------------------------------------
 
@@ -1046,9 +1101,11 @@
          WRITE(numout,*) '         Namelist nam_zgr_sco'
          WRITE(numout,*) '            sigma-stretching coeffs '
-         WRITE(numout,*) '            maximum depth of s-bottom surface (>0)    sbot_max = ' ,sbot_max
-         WRITE(numout,*) '            minimum depth of s-bottom surface (>0)    sbot_min = ' ,sbot_min
-         WRITE(numout,*) '            surface control parameter (0<=theta<=20)  theta    = ', theta
-         WRITE(numout,*) '            bottom control parameter  (0<=thetb<= 1)  thetb    = ', thetb
-         WRITE(numout,*) '            maximum cut-off r-value allowed           r_max    = ' , r_max
+         WRITE(numout,*) '            maximum depth of s-bottom surface (>0)    sbot_max   = ' ,sbot_max
+         WRITE(numout,*) '            minimum depth of s-bottom surface (>0)    sbot_min   = ' ,sbot_min
+         WRITE(numout,*) '            surface control parameter (0<=theta<=20)  theta      = ', theta
+         WRITE(numout,*) '            bottom control parameter  (0<=thetb<= 1)  thetb      = ', thetb
+         WRITE(numout,*) '            maximum cut-off r-value allowed           r_max      = ' , r_max
+         WRITE(numout,*) '            Critical depth                            hc         = ', hc
+         WRITE(numout,*) '            Hybrid s-sigma-coordinate                 ln_s_sigma = ', ln_s_sigma
       ENDIF
 
@@ -1057,5 +1114,6 @@
       hifv(:,:) = sbot_min
       hiff(:,:) = sbot_min
-      !                                        ! set maximum ocean depth 
+
+      !                                        ! set maximum ocean depth
       bathy(:,:) = MIN( sbot_max, bathy(:,:) )
 
@@ -1063,7 +1121,7 @@
       !                                        ! Define the envelop bathymetry   (hbatt)
       !                                        ! =============================
-      ! Smooth the bathymetry (if required) 
+      ! Smooth the bathymetry (if required)
       scosrf(:,:) = 0.e0              ! ocean surface depth (here zero: no under ice-shelf sea)
-      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth 
+      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
       !
       ! use r-value to create hybrid coordinates
@@ -1099,8 +1157,8 @@
          DO jj = 1, nlcj
             DO ji = 1, nlci
-                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )    
-            END DO
-         END DO
-         ! 
+                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(:,:) ) )
          !
@@ -1124,5 +1182,5 @@
             END DO
          END DO
-         ! 
+         !
          DO jj = 1, nlcj
             DO ji = 1, nlci
@@ -1237,19 +1295,95 @@
       !
       ! non-dimensional "sigma" for model level depth at w- and t-levels
-      DO jk = 1, jpk
-        gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
-        gsigt(jk) = -fssig( REAL(jk,wp)        )
-      END DO
-      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
-      !
-      ! Coefficients for vertical scale factors at w-, t- levels
+
+      IF ( ln_s_sigma ) THEN  !Song and Haidvogel style stretched sigma for depths below hc, with uniform sigma in shallower waters
+
+         DO ji=1,jpi
+           DO jj=1,jpj
+
+             IF (hbatt(ji,jj).GT.hc) THEN !deep water, stretched sigma
+               DO jk = 1, jpk
+                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, bb )
+                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , bb )
+               END DO
+             ELSE ! shallow water, uniform sigma
+               DO jk = 1, jpk
+                 gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
+                 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
+               END DO
+             ENDIF
+             IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
+
+
+             DO jk = 1, jpkm1
+                esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
+                esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
+             END DO
+             esigw3(ji,jj,1) = esigw3(ji,jj,2)
+             esigt3(ji,jj,jpk) = esigt3(ji,jj,jpkm1)
+
+             ! Coefficients for vertical depth as the sum of e3w scale factors
+             gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
+             DO jk = 2, jpk
+                gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
+             END DO
+
+             DO jk = 1, jpk
+                zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp)
+                zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp)
+                gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsigt3(ji,jj,jk)+hc*zcoeft)
+                gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsigw3(ji,jj,jk)+hc*zcoefw)
+                gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsi3w3(ji,jj,jk)+hc*zcoefw)
+             END DO
+
+           ENDDO   ! for all jj's
+         ENDDO    ! for all ji's
+
+
+         DO ji=1,jpi
+           DO jj=1,jpj
+
+             DO jk = 1, jpk
+                esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / &
+                                   ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+                esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / &
+                                   ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+                esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  &
+                                     hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / &
+                                   ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+                esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / &
+                                   ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+                esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / &
+                                   ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+
+                e3t(ji,jj,jk)=((hbatt(ji,jj)-hc)*esigt3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                e3u(ji,jj,jk)=((hbatu(ji,jj)-hc)*esigtu3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                e3v(ji,jj,jk)=((hbatv(ji,jj)-hc)*esigtv3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                e3f(ji,jj,jk)=((hbatf(ji,jj)-hc)*esigtf3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                !
+                e3w (ji,jj,jk)=((hbatt(ji,jj)-hc)*esigw3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                e3uw(ji,jj,jk)=((hbatu(ji,jj)-hc)*esigwu3(ji,jj,jk) + hc/FLOAT(jpkm1))
+                e3vw(ji,jj,jk)=((hbatv(ji,jj)-hc)*esigwv3(ji,jj,jk) + hc/FLOAT(jpkm1))
+             END DO
+
+           ENDDO
+         ENDDO
+
+      ELSE   ! not ln_s_sigma
+
+         DO jk = 1, jpk
+           gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
+           gsigt(jk) = -fssig( REAL(jk,wp)     )
+         END DO
+         IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
+         !
+     ! Coefficients for vertical scale factors at w-, t- levels
 !!gm bug :  define it from analytical function, not like juste bellow....
 !!gm        or betteroffer the 2 possibilities....
-      DO jk = 1, jpkm1
-         esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
-         esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
-      END DO
-      esigw( 1 ) = esigw(  2  )
-      esigt(jpk) = esigt(jpkm1)
+         DO jk = 1, jpkm1
+            esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
+            esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
+         END DO
+         esigw( 1 ) = esigw(  2  )
+         esigt(jpk) = esigt(jpkm1)
 
 !!gm  original form
@@ -1259,33 +1393,37 @@
 !!org END DO
 !!gm
-      !
-      ! Coefficients for vertical depth as the sum of e3w scale factors
-      gsi3w(1) = 0.5 * esigw(1)
-      DO jk = 2, jpk
-         gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
-      END DO
+         !
+         ! Coefficients for vertical depth as the sum of e3w scale factors
+         gsi3w(1) = 0.5 * esigw(1)
+         DO jk = 2, jpk
+            gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
+         END DO
 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
-      DO jk = 1, jpk
-         zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
-         zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 
-         gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)
-         gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)
-         gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw)
-      END DO
+         DO jk = 1, jpk
+            zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
+            zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
+            gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)
+            gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)
+            gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw)
+         END DO
 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            DO jk = 1, jpk
-              e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
-              e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
-              e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
-              e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
-              !
-              e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
-              e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
-              e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
-            END DO
-         END DO
-      END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 1, jpk
+                 e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
+                 e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
+                 e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
+                 e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
+                 !
+                 e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
+                 e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
+                 e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
+               END DO
+            END DO
+         END DO
+
+      ENDIF ! ln_s_sigma
+
+
       !
       ! HYBRID : 
