Index: branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 2864)
+++ branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 2908)
@@ -21,5 +21,5 @@
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
-  !!   dia_dct      :  compute the transport trough a sec.
+  !!   dia_dct      :  compute the transport through a sec.
   !!   dia_dct_init :  read namelist.
   !!   readsec      :  read sections description and pathway
@@ -78,22 +78,22 @@
 
   TYPE SECTION
-     CHARACTER(len=60)                            :: name           ! name of the sec
-     LOGICAL                                      :: llstrpond      ! true if you want the computation of salinity and
-                                                                    ! temperature balanced by the transport
-     LOGICAL                                      :: ll_ice_section ! icesurf and icevol computation
-     LOGICAL                                      :: ll_date_line   ! = T if the section crosses the date-line
-     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec       ! longitude and latitude of the extremities of the sec
-     INTEGER                                      :: nb_class       ! number of boundaries for density classes
-     INTEGER, DIMENSION(nb_point_max)             :: direction      ! vector direction of the point in the section
-     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname      ! caracteristics of the class
-     REAL(wp), DIMENSION(nb_class_max)            :: zsigi        ,&! in-situ   density classes    (99 if you don't want)
-                                                     zsigp        ,&! potential density classes    (99 if you don't want)
-                                                     zsal         ,&! salinity classes   (99 if you don't want)
-                                                     ztem         ,&! temperature classes(99 if you don't want)
-                                                     zlay           ! level classes      (99 if you don't want)
+     CHARACTER(len=60)                            :: name              ! name of the sec
+     LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
+                                                                       ! heat transports
+     LOGICAL                                      :: ll_ice_section    ! icesurf and icevol computation
+     LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
+     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
+     INTEGER                                      :: nb_class          ! number of boundaries for density classes
+     INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
+     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! caracteristics of the class
+     REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
+                                                     zsigp           ,&! potential density classes    (99 if you don't want)
+                                                     zsal            ,&! salinity classes   (99 if you don't want)
+                                                     ztem            ,&! temperature classes(99 if you don't want)
+                                                     zlay              ! level classes      (99 if you don't want)
      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
-     REAL(wp)                                         :: slopeSection  ! coeff directeur de la section
-     INTEGER                                          :: nb_point              ! nb de points de la section
-     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint             ! list of point in the section
+     REAL(wp)                                         :: slopeSection  ! section's slopesection 
+     INTEGER                                          :: nb_point      ! section's number of points
+     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! section's list point
   END TYPE SECTION
 
@@ -161,4 +161,11 @@
      LOGICAL             :: lldebug =.FALSE.  !debug a section  
      CHARACTER(len=160)  :: clfileout         !fileout name
+
+     
+     INTEGER , DIMENSION(1):: ish                                       ! tmp array for mpp_sum
+     INTEGER , DIMENSION(3):: ish2                                      !   "
+     REAL(wp), DIMENSION(nb_sec_max*nb_type_class*nb_class_max):: zwork !   "  
+     REAL(wp), DIMENSION(nb_sec_max,nb_type_class,nb_class_max):: zsum  !   "
+
      !!---------------------------------------------------------------------    
 
@@ -183,10 +190,25 @@
            !Compute transport through section  
            CALL transport(secs(jsec),lldebug) 
+
+        ENDDO
              
-           IF( MOD(kt,nn_dctwri)==0 )THEN
-
-              IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: write at kt = ",kt         
+        IF( MOD(kt,nn_dctwri)==0 )THEN
+
+           IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: write at kt = ",kt         
   
-              !Write the transport
+           !Sum on all procs 
+           IF( lk_mpp )THEN
+              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
+              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
+              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
+              zwork(:)= RESHAPE(zsum(:,:,:), ish )
+              CALL mpp_sum(zwork, ish(1))
+              zsum(:,:,:)= RESHAPE(zwork,ish2)
+              DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
+           ENDIF
+
+           !Write the transport
+           DO jsec=1,nb_sec
+
               IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec))
             
@@ -194,6 +216,8 @@
               secs(jsec)%transport(:,:)=0.  
 
-           ENDIF 
-        ENDDO
+           ENDDO
+
+        ENDIF 
+
      ENDIF
 
@@ -310,9 +334,4 @@
               ENDDO                  
            ENDIF
-           !?IF( jsec==nn_secdebug .OR. nn_secdebug==-1  )THEN
-           !  DO jpt=1,iptglo
-           !      WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt)
-           !   ENDDO                  
-           !ENDIF
  
            !Now each proc selects only points that are in its domain:
@@ -335,5 +354,4 @@
                  secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
                  secs(jsec)%direction(iptloc) = directemp(jpt)               ! store local direction
-                 !WRITE(narea+200,*)'         # I J : ',iiloc,ijloc
               ENDIF
 
@@ -356,7 +374,6 @@
 
            !remove redundant points between processors
-           !look NEMO documentation for more explanations
-!           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
-           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1)  ) lldebug = .TRUE.
+           !------------------------------------------
+           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
            IF( iptloc .NE. 0 )THEN
               CALL removepoints(secs(jsec),'I','top_list',lldebug)
@@ -409,4 +426,5 @@
                                  !where points could be redundant
                 isgn          ,& !way of course in listpoint
+                ipoint          ,& !way of course in listpoint
                 istart,iend      !first and last points selected in listpoint
      INTEGER :: jpoint   =0      !loop on list points
@@ -443,8 +461,9 @@
 
      jpoint=iextr+isgn
-     DO WHILE ( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point .AND. &
-             icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest)
-             jpoint=jpoint+isgn
+     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point  .AND. &
+        icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )
+        jpoint=jpoint+isgn
      ENDDO
+
      IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
      ELSE                        ; istart=1        ; iend=jpoint+1
@@ -489,5 +508,4 @@
      !! * Local variables
      INTEGER             :: jk,jseg,jclass,   &!loop on level/segment/classes 
-                            iptegrid,         &!=0 if section is constant along i/j 
                             isgnu  , isgnv     !
      INTEGER :: ii, ij ! local integer
@@ -506,8 +524,5 @@
 
      TYPE(POINT_SECTION) :: k
-     INTEGER,DIMENSION(1):: ish
-     INTEGER,DIMENSION(2):: ish2
      REAL(wp),DIMENSION(nb_type_class,nb_class_max)::zsum
-     REAL(wp),DIMENSION(nb_type_class*nb_class_max)::zwork  ! temporary vector for mpp_sum
      !!--------------------------------------------------------
 
@@ -520,13 +535,4 @@
      zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp
      zice_vol_pos  = 0._wp ; zice_vol_neg  = 0._wp
-
-     !compute iptegrid: iptegrid=0 if section is constant along i or j
-     iptegrid = 0
-     jseg     = 1
-     DO WHILE( (jseg .LT. MAX(sec%nb_point-1,0)) .AND. iptegrid .NE. 1 )
-        IF(sec%direction(jseg) .NE. sec%direction(1)) iptegrid = 1    
-        jseg = jseg+1
-     ENDDO 
-     IF(lk_mpp) CALL mpp_sum(iptegrid)
 
      !---------------------------!
@@ -623,5 +629,4 @@
                  !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 
                  !----------------------------------------------!
-           IF(ld_debug)write(narea+200,*)k ,jk
  
                  IF ( (    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.    &
@@ -663,5 +668,4 @@
                     ENDIF
 #endif
-           !IF(ld_debug)write(narea+200,*)k ,jk,zTnorm             
                     !COMPUTE TRANSPORT 
                     !zTnorm=transport through one cell for one class
@@ -771,14 +775,4 @@
      ENDIF   !end of nb_inmesh=0 case
 
-     !------------------!
-     !SUM ON ALL PROCS  !
-     !------------------!
-     IF( lk_mpp )THEN
-         ish(1)  =  nb_type_class*nb_class_max ;  ish2(1)=nb_type_class ; ish2(2)=nb_class_max
-         zwork(:)= RESHAPE(zsum(:,:), ish )
-         CALL mpp_sum(zwork, ish(1))
-         zsum(:,:)= RESHAPE(zwork,ish2)
-     ENDIF
-     
      !-------------------------------|
      !FINISH COMPUTING TRANSPORTS    |
@@ -838,27 +832,30 @@
      !! Method:
      !!        1. Write volume transports in "volume_transport"
-     !!           Unit: Sv : Surface * Velocity / 1.e6 
+     !!           Unit: Sv : area * Velocity / 1.e6 
      !! 
      !!        2. Write heat transports in "heat_transport"
-     !!           Unit: Peta W : Surface * Velocity * T * rhau * Cp / 1.e15
+     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
      !! 
      !!        3. Write salt transports in "salt_transport"
-     !!           Unit: Mega m^3 / s : Surface * Velocity * S / 1.e6
+     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
      !!
      !!------------------------------------------------------------- 
      !!arguments
-     INTEGER, INTENT(IN)          :: kt
-     TYPE(SECTION), INTENT(INOUT) :: sec
-     INTEGER ,INTENT(IN)          :: jsec
+     INTEGER, INTENT(IN)          :: kt         ! time-step
+     TYPE(SECTION), INTENT(INOUT) :: sec        ! section to write   
+     INTEGER ,INTENT(IN)          :: jsec       ! section number
 
      !!local declarations
      REAL(wp) ,DIMENSION(nb_type_class):: zsumclass
-     INTEGER               :: jcl,ji        !Dummy loop
-     CHARACTER(len=2)      :: classe
-     REAL(wp)              :: zbnd1,zbnd2
+     INTEGER               :: jcl,ji            ! Dummy loop
+     CHARACTER(len=2)      :: classe            ! Classname 
+     REAL(wp)              :: zbnd1,zbnd2       ! Class bounds
+     REAL(wp)              :: zslope            ! section's slope coeff
      !!------------------------------------------------------------- 
       
      zsumclass(:)=0._wp
-               
+     zslope = sec%slopeSection       
+
+ 
      DO jcl=1,MAX(1,sec%nb_class-1)
 
@@ -900,6 +897,6 @@
         ENDIF
         !temperature classes transports
-        IF( ( sec%ztem(jcl+1) .NE. 99.     ) .AND. &
-            ( sec%ztem(jcl+1) .NE. 99.     )       ) THEN
+        IF( ( sec%ztem(jcl+1) .NE. 99._wp     ) .AND. &
+            ( sec%ztem(jcl+1) .NE. 99._wp     )       ) THEN
            classe = 'T       '
            zbnd1 = sec%ztem(jcl)
@@ -908,5 +905,5 @@
                   
         !write volume transport per class
-        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, &
+        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope, &
                               jcl,classe,zbnd1,zbnd2,&
                               sec%transport(1,jcl),sec%transport(2,jcl), &
@@ -916,10 +913,10 @@
 
            !write heat transport per class:
-           WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  &
+           WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,zslope,  &
                               jcl,classe,zbnd1,zbnd2,&
                               sec%transport(7,jcl)*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, &
                               ( sec%transport(7,jcl)+sec%transport(8,jcl) )*1000._wp*rcp/1.e15
            !write salt transport per class
-           WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  &
+           WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,zslope,  &
                               jcl,classe,zbnd1,zbnd2,&
                               sec%transport(9,jcl)*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,&
@@ -934,5 +931,5 @@
 
      !write total volume transport
-     WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, &
+     WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope, &
                            jcl,"total",zbnd1,zbnd2,&
                            zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2)
@@ -941,11 +938,13 @@
 
         !write total heat transport
-        WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
+        WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,zslope, &
                            jcl,"total",zbnd1,zbnd2,&
-                           zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,(zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15
+                           zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,&
+                           (zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15
         !write total salt transport
-        WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
+        WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,zslope, &
                            jcl,"total",zbnd1,zbnd2,&
-                           zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,(zsumclass(9)+zsumclass(10))*1000._wp/1.e9
+                           zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,&
+                           (zsumclass(9)+zsumclass(10))*1000._wp/1.e9
      ENDIF
 
@@ -953,10 +952,10 @@
      IF ( sec%ll_ice_section) THEN
         !write total ice volume transport
-        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
+        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope,&
                               jcl,"ice_vol",zbnd1,zbnd2,&
                               sec%transport(9,1),sec%transport(10,1),&
                               sec%transport(9,1)+sec%transport(10,1)
         !write total ice surface transport
-        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
+        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope,&
                               jcl,"ice_surf",zbnd1,zbnd2,&
                               sec%transport(11,1),sec%transport(12,1), &
@@ -972,38 +971,72 @@
   !!----------------------------------------------------------------------
   !!
-  !!   Purpose: compute Temperature/Salinity/density on U-point or V-point
+  !!   Purpose: compute Temperature/Salinity/density at U-point or V-point
+  !!   --------
   !!
-  !!    |    I       |    I+1     |
-  !!    |            |            |
-  !!  --------------------------------     1.Compute Tbis:
-  !!    |            |            |          interpolate T(I,J,K) et T(I,J,K+1)
-  !!    |            |            |
-  !!    |            |            |        2. 1.Compute Temp/Salt/Rho à U point: 
-  !! K  |    T       |            |          interpolate Tbis et T(I+U,J,K+1)  
-  !!    |            |            |
-  !!    |            |            |
-  !!    |            |            |
-  !!  --------------------------------
-  !!    |            |            |
-  !!    |            |            |
-  !!    |            |            |
-  !!K+1 |    Tbis    U     T      |
-  !!    |            |            |
-  !!    |    T       |            |
-  !!    |            |------------|
-  !!    |            | partials   |
-  !!    |            |  steps     |
-  !!  --------------------------------
+  !!   Method:
+  !!   ------
+  !!
+  !!   ====> full step and partial step
+  !! 
+  !!
+  !!    |    I          |    I+1           |    Z=Temperature/Salinity/density at U-poinT
+  !!    |               |                  |
+  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis
+  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
+  !!    |               |                  |       zbis = 
+  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
+  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ] 
+  !!    |               |                  | 
+  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
+  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K)  
+  !!    |     .         |                  |
+  !!    |     .         |                  |       Z= ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1) 
+  !!    |     .         |                  |
+  !!  ------------------------------------------
+  !!    |     .         |                  |
+  !!    |     .         |                  |
+  !!    |     .         |                  |
+  !!K   |    zbis...... U ..ptab(I+1,J,K)  |
+  !!    |     .         |                  |
+  !!    | ptab(I,J,K)   |                  |
+  !!    |               |------------------|
+  !!    |               | partials         |
+  !!    |               |  steps           |
+  !!  -------------------------------------------
+  !!    <----zet1------><----zet2--------->
+  !!
+  !!
+  !!   ====>  s-coordinate
+  !!     
+  !!    |                |                  |   1. Compute distance between T1 and U points
+  !!    |                |                  |      Compute distance between T2 and U points
+  !!    |                | ptab(I+1,J,K)    | 
+  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point 
+  !!    |                |      ^           |    
+  !!    |                |      | zdep2     |    
+  !!    |                |      |           |    
+  !!    |       ^        U      v           |
+  !!    |       |        |                  |
+  !!    |       | zdep1  |                  |    
+  !!    |       v        |                  |
+  !!    |      T1        |                  |
+  !!    | ptab(I,J,K)    |                  | 
+  !!    |                |                  | 
+  !!    |                |                  | 
+  !!
+  !!    <----zet1--------><----zet2--------->
   !!
   !!----------------------------------------------------------------------
   !*arguments
-  INTEGER, INTENT(IN)                          :: ki, kj, kk
-  CHARACTER(len=1), INTENT(IN)                 :: cd_point
-  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab
-  REAL(wp)                                     :: interp
+  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
+  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
+  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
+  REAL(wp)                                     :: interp       ! interpolated variable 
 
   !*local declations
-  INTEGER :: ii1, ij1, ii2, ij2
-  REAL(wp):: ze3w, zfse3, zmax1, zmax2, zbis
+  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
+  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
+  REAL(wp):: zet1, zet2                                        ! weight for interpolation 
+  REAL(wp):: zdep1,zdep2                                       ! differences of depth
   !!----------------------------------------------------------------------
 
@@ -1011,44 +1044,70 @@
      ii1 = ki    ; ij1 = kj 
      ii2 = ki+1  ; ij2 = kj 
+
+     zet1=e1t(ii1,ij1)
+     zet2=e1t(ii2,ij2)
+
   ELSE ! cd_point=='V' 
      ii1 = ki    ; ij1 = kj 
      ii2 = ki    ; ij2 = kj+1  
+
+     zet1=e2t(ii1,ij1)
+     zet2=e2t(ii2,ij2)
+
   ENDIF
 
+  IF( ln_sco )THEN   ! s-coordinate case
+
+     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2 
+     zdep1 = fsdept(ii1,ij1,kk) - zdepu
+     zdep2 = fsdept(ii2,ij2,kk) - zdepu
+
+     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
+     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
+  
+     ! result
+     interp = umask(ii1,ij1,kk) * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )   
+
+
+  ELSE       ! full step or partial step case 
+
 #if defined key_vvl
 
-  ze3w  = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk)
-
-  zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) )
-  zmax1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk)
+     ze3w  = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk)
+
+     zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) )
+     zwgt1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk)
   
-  zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) )
-  zmax2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk)
+     zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) )
+     zwgt2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk)
 
 #else
 
-  ze3w  = fse3w(ii2,ij2,kk)-fse3w(ii1,ij1,kk) 
-  zmax1 =  ze3w / fse3w(ii2,ij2,kk)
-  zmax2 = -ze3w / fse3w(ii1,ij1,kk)
+     ze3t  = fse3t(ii2,ij2,kk)-fse3t(ii1,ij1,kk) 
+     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
+     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
 
 #endif
 
-  IF(kk .NE. 1)THEN
-
-     IF( ze3w >= 0. )THEN 
-        !zbis
-        zbis = ptab(ii2,ij2,kk) + zmax1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
-        ! result
-        interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + zbis )
+     IF(kk .NE. 1)THEN
+
+        IF( ze3t >= 0. )THEN 
+           !zbis
+           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
+           ! result
+            interp = umask(ii1,ij1,kk) * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
+        ELSE
+           !zbis
+           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
+           ! result
+           interp = umask(ii1,ij1,kk) * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
+        ENDIF    
+
      ELSE
-        !zbis
-        zbis = ptab(ii1,ij1,kk) + zmax2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
-        ! result
-        interp = umask(ii1,ij1,kk) * 0.5*( zbis +  ptab(ii2,ij2,kk) )
-     ENDIF    
-
-  ELSE
-     interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + ptab(ii2,ij2,kk) )
+        interp = umask(ii1,ij1,kk) * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
+     ENDIF
+
   ENDIF
+
 
   END FUNCTION interp
