Index: /branches/2012/dev_r3521_INGV7_muscl/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
===================================================================
--- /branches/2012/dev_r3521_INGV7_muscl/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 3526)
+++ /branches/2012/dev_r3521_INGV7_muscl/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 3527)
@@ -8,4 +8,5 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
+   !!            3.4  !  2012-06  (P. Oddo) include the upstream where needed
    !!----------------------------------------------------------------------
 
@@ -18,7 +19,9 @@
    USE trdmod_oce      ! tracers trends 
    USE trdtra      ! tracers trends 
+   USE eosbn2          ! equation of state 
    USE in_out_manager  ! I/O manager
    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE trabbl          ! tracers: bottom boundary layer
+   USE sbcrnf          ! river runoffs
    USE lib_mpp         ! distribued memory computing
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
@@ -27,4 +30,6 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE eosbn2          ! equation of state
+   USE sbcrnf          ! river runoffs
 
    IMPLICIT NONE
@@ -35,4 +40,6 @@
    LOGICAL  :: l_trd       ! flag to compute trends
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
+   !                                                             !  and in closed seas (orca 2 and 4 configurations)
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -74,8 +81,13 @@
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      INTEGER  ::   ierr
       REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars
       REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      -
       REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy
+      INTEGER  ::   ierr
+      REAL(wp) ::   zice                             ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztfreez
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zind
       !!----------------------------------------------------------------------
       !
@@ -83,4 +95,6 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy )
+      CALL wrk_alloc( jpi, jpj, ztfreez )
+      CALL wrk_alloc( jpi, jpj, jpk, zind )
       !
 
@@ -89,4 +103,13 @@
          IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         !
+         !
+         IF (.not. ALLOCATED(upsmsk))THEN
+             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
+             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate array')
+         ENDIF
+         !
+         upsmsk(:,:) = 0._wp                             ! not upstream by default
          !
          l_trd = .FALSE.
@@ -94,4 +117,24 @@
       ENDIF
 
+      !
+      ! Upstream / centered scheme indicator
+      ! ------------------------------------
+      ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
+               IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1_wp ) THEN   ;   zice = 1.e0
+               ELSE                                                         ;   zice = 0.e0
+               ENDIF
+               zind(ji,jj,jk) = MAX (   &
+                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
+                  upsmsk(ji,jj)               ,      &  ! some of some straits
+                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
+                  &                  ) * tmask(ji,jj,jk)
+               zind(ji,jj,jk) = 1 - zind(ji,jj,jk)
+            END DO
+         END DO
+      END DO
       !                                                     ! ===========
       DO jn = 1, kjpt                                       ! tracer loop
@@ -148,6 +191,6 @@
                   zalpha = 0.5 - z0u
                   zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
-                  zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk)
-                  zzwy = ptb(ji  ,jj,jk,jn) + zu * zslpx(ji  ,jj,jk)
+                  zzwx = ptb(ji+1,jj,jk,jn) + zind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk))
+                  zzwy = ptb(ji  ,jj,jk,jn) + zind(ji,jj,jk) * (zu * zslpx(ji  ,jj,jk))
                   zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                   !
@@ -155,6 +198,6 @@
                   zalpha = 0.5 - z0v
                   zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
-                  zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk)
-                  zzwy = ptb(ji,jj  ,jk,jn) + zv * zslpy(ji,jj  ,jk) 
+                  zzwx = ptb(ji,jj+1,jk,jn) + zind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk))
+                  zzwy = ptb(ji,jj  ,jk,jn) + zind(ji,jj,jk) * (zv * zslpy(ji,jj  ,jk))
                   zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                END DO
@@ -230,6 +273,6 @@
                   zalpha = 0.5 + z0w
                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 
-                  zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1)
-                  zzwy = ptb(ji,jj,jk  ,jn) + zw * zslpx(ji,jj,jk  )
+                  zzwx = ptb(ji,jj,jk+1,jn) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1))
+                  zzwy = ptb(ji,jj,jk  ,jn) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk  ))
                   zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                END DO 
@@ -255,4 +298,6 @@
       !
       CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy )
+      CALL wrk_dealloc( jpi, jpj, ztfreez )
+      CALL wrk_dealloc( jpi, jpj, jpk, zind )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl')
