Index: branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
===================================================================
--- branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 5947)
+++ branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 5948)
@@ -16,6 +16,6 @@
    USE trc_oce        ! share passive tracers/Ocean variables
    USE trd_oce        ! trends: ocean variables
+   USE traadv_fct      ! acces to routine interp_4th_cpt 
    USE trdtra         ! trends manager: tracers 
-   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient
    USE diaptr         ! poleward transport diagnostics
    !
@@ -38,5 +38,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -44,6 +44,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
-      &                                       ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn,          &
+      &                                                ptb, ptn, pta, kjpt, kn_ubs_v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_ubs  ***
@@ -52,5 +52,5 @@
       !!      and add it to the general trend of passive tracer equations.
       !!
-      !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order
+      !! ** Method  :   The 3rd order Upstream Biased Scheme (UBS) is based on an
       !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
       !!      It is only used in the horizontal direction.
@@ -61,14 +61,17 @@
       !!      where zltu is the second derivative of the before temperature field:
       !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
-      !!      This results in a dissipatively dominant (i.e. hyper-diffusive) 
+      !!        This results in a dissipatively dominant (i.e. hyper-diffusive) 
       !!      truncation error. The overall performance of the advection scheme 
       !!      is similar to that reported in (Farrow and Stevens, 1995). 
-      !!      For stability reasons, the first term of the fluxes which corresponds
+      !!        For stability reasons, the first term of the fluxes which corresponds
       !!      to a second order centered scheme is evaluated using the now velocity 
       !!      (centered in time) while the second term which is the diffusive part 
       !!      of the scheme, is evaluated using the before velocity (forward in time). 
       !!      Note that UBS is not positive. Do not use it on passive tracers.
-      !!                On the vertical, the advection is evaluated using a TVD scheme,
-      !!      as the UBS have been found to be too diffusive.
+      !!                On the vertical, the advection is evaluated using a FCT scheme,
+      !!      as the UBS have been found to be too diffusive. 
+!!gm  !!                kn_ubs_v argument (not coded for the moment)
+      !!      controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2) 
+      !!      or on a 4th order compact scheme (kn_ubs_v=4).
       !!
       !! ** Action : - update (pta) with the now advective tracer trends
@@ -81,4 +84,5 @@
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
+      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
       REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components
@@ -95,5 +99,5 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw )
+      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zltu, zltv, zti, ztw )
       !
       IF( kt == kit000 )  THEN
@@ -106,25 +110,24 @@
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
       !
+      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp     ! Bottom value : set to zero one for all
+      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp
+      IF( lk_vvl )   ztw(:,:, 1 ) = 0._wp                   ! surface value: set to zero only in vvl case
+      !
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         ! 1. Bottom value : flux set to zero
-         ! ----------------------------------
-         zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0
          !                                              
-         DO jk = 1, jpkm1                                 ! Horizontal slab
-            !                                   
-            !  Laplacian
-            DO jj = 1, jpjm1            ! First derivative (gradient)
+         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
+            DO jj = 1, jpjm1              ! First derivative (masked gradient)
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
-                  zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
+                  zeeu = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
+                  zeev = e1_e2v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)
                   ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
                   ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
                END DO
             END DO
-            DO jj = 2, jpjm1            ! Second derivative (divergence)
+            DO jj = 2, jpjm1              ! Second derivative (divergence)
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
+                  zcoef = 1._wp / ( 6._wp * fse3t(ji,jj,jk) )
                   zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
                   zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
@@ -132,30 +135,27 @@
             END DO
             !                                    
-         END DO                                           ! End of slab         
+         END DO         
          CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
-
          !    
-         !  Horizontal advective fluxes               
-         DO jk = 1, jpkm1                                 ! Horizontal slab
+         DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS)
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  ! upstream transport (x2)
-                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
+                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2)
                   zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
                   zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
                   zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
-                  ! 2nd order centered advective fluxes (x2)
+                  !                                                  ! 2nd order centered advective fluxes (x2)
                   zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
                   zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
-                  ! UBS advective fluxes
+                  !                                                  ! UBS advective fluxes
                   ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
                   ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
                END DO
             END DO
-         END DO                                           ! End of slab         
-
-         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends
-
-         DO jk = 1, jpkm1                 ! Horizontal advective trends
+         END DO         
+         !
+         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update
+         !
+         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -166,9 +166,8 @@
             END DO
             !                                             
-         END DO                                           !   End of slab
-
-         ! Horizontal trend used in tra_adv_ztvd subroutine
-         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
-
+         END DO
+         !
+         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
+         !                                            ! and/or in trend diagnostic (l_trd=T) 
          !                
          IF( l_trd ) THEN                  ! trend diagnostics
@@ -181,76 +180,89 @@
             IF( jn == jp_sal )  str_adv(:) = ptr_sj( ztv(:,:,:) )
          ENDIF
-         
-         ! TVD scheme for the vertical direction  
-         ! ----------------------
-         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
-
-         !  Bottom value : flux set to zero
-         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
-
-         ! Surface value
-         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
-         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface 
-         ENDIF
-         !  upstream advection with initial mass fluxes & intermediate update
-         ! -------------------------------------------------------------------
-         ! Interior value
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
-                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
-                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
-               END DO
-            END DO
-         END DO 
-         ! update and guess with monotonic sheme
-         DO jk = 1, jpkm1
-            z2dtt = p2dt(jk)
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
-                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
-                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         !
-         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
-
-         !  antidiffusive flux : high order minus low order
-         ztw(:,:,1) = 0.e0       ! Surface value
-         DO jk = 2, jpkm1        ! Interior value
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         !
-         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
-
-         !  final trend with corrected fluxes
-         DO jk = 1, jpkm1
+         !
+         !                       !== vertical advective trend  ==!
+         !
+         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme
+         !
+         CASE(  2  )                   ! 2nd order FCT 
+            !         
+            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
+            !
+            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
+            DO jk = 2, jpkm1                 ! Interior value (w-masked)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
+                     zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
+                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO 
+            IF(.NOT.lk_vvl ) THEN            ! top ocean value (only in linear free surface as ztw has been w-masked)
+               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                     END DO
+                  END DO   
+               ELSE                                ! no cavities: only at the ocean surface
+                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+               ENDIF
+            ENDIF
+            !
+            DO jk = 1, jpkm1           !* trend and after field with monotonic scheme
+               z2dtt = p2dt(jk)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
+                     pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
+                     zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
+            !
+            !                          !*  anti-diffusive flux : high order minus low order
+            DO jk = 2, jpkm1        ! Interior value  (w-masked)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
+                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0
+            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = 0._wp      ! only ocean surface as interior zwz values have been w-masked
+            !
+            CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
+            !
+         CASE(  4  )                               ! 4th order COMPACT
+            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point
+            DO jk = 2, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            IF(.NOT.lk_vvl )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
+            !
+         END SELECT
+         !
+         DO jk = 1, jpkm1        !  final trend with corrected fluxes
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.   
-                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ! k- vertical advective trends  
-                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
-                  ! added to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
-               END DO
-            END DO
-         END DO
-
-         !  Save the final vertical advective trends
-         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         !
+         IF( l_trd )  THEN       ! vertical advective trend diagnostics
             DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
-                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
+                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          &
+                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   &
+                        &                              / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                   END DO
                END DO
@@ -261,5 +273,5 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw )
+      CALL wrk_dealloc( jpi,jpj,jpk,   ztu, ztv, zltu, zltv, zti, ztw )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs')
@@ -294,17 +306,17 @@
       IF( nn_timing == 1 )  CALL timing_start('nonosc_z')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo )
+      CALL wrk_alloc( jpi,jpj,jpk,   zbetup, zbetdo )
       !
       zbig  = 1.e+40_wp
       zrtrn = 1.e-15_wp
       zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
-
+      !
       ! Search local extrema
       ! --------------------
-      ! large negative value (-zbig) inside land
+      !                    ! large negative value (-zbig) inside land
       pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
       paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
-      ! search maximum in neighbourhood
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1     ! search maximum in neighbourhood
          ikm1 = MAX(jk-1,1)
          DO jj = 2, jpjm1
@@ -316,9 +328,9 @@
          END DO
       END DO
-      ! large positive value (+zbig) inside land
+      !                    ! large positive value (+zbig) inside land
       pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
       paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
-      ! search minimum in neighbourhood
-      DO jk = 1, jpkm1
+      !
+      DO jk = 1, jpkm1     ! search minimum in neighbourhood
          ikm1 = MAX(jk-1,1)
          DO jj = 2, jpjm1
@@ -330,13 +342,10 @@
          END DO
       END DO
-
-      ! restore masked values to zero
+      !                    ! restore masked values to zero
       pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
       paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
-
-
-      ! 2. Positive and negative part of fluxes and beta terms
-      ! ------------------------------------------------------
-
+      !
+      ! Positive and negative part of fluxes and beta terms
+      ! ---------------------------------------------------
       DO jk = 1, jpkm1
          z2dtt = p2dt(jk)
@@ -347,5 +356,5 @@
                zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
                ! up & down beta terms
-               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
+               zbt = e1e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
                zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
                zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
@@ -353,4 +362,5 @@
          END DO
       END DO
+      !
       ! monotonic flux in the k direction, i.e. pcc
       ! -------------------------------------------
@@ -366,5 +376,5 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zbetup, zbetdo )
       !
       IF( nn_timing == 1 )  CALL timing_stop('nonosc_z')
