Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/daymod.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/daymod.F90	(revision 6879)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/daymod.F90	(revision 6879)
@@ -0,0 +1,418 @@
+MODULE daymod
+   !!======================================================================
+   !!                       ***  MODULE  daymod  ***
+   !! Ocean        :  calendar
+   !!=====================================================================
+   !! History :  OPA  ! 1994-09  (M. Pontaud M. Imbard)  Original code
+   !!                 ! 1997-03  (O. Marti)
+   !!                 ! 1997-05  (G. Madec)
+   !!                 ! 1997-08  (M. Imbard)
+   !!   NEMO     1.0  ! 2003-09  (G. Madec)  F90 + nyear, nmonth, nday
+   !!                 ! 2004-01  (A.M. Treguier) new calculation based on adatrj
+   !!                 ! 2006-08  (G. Madec)  surface module major update
+   !!                 ! 2015-11  (D. Lea) Allow non-zero initial time of day
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   day        : calendar
+   !!
+   !!           -------------------------------
+   !!           ----------- WARNING -----------
+   !!
+   !!   we suppose that the time step is deviding the number of second of in a day
+   !!             ---> MOD( rday, rdt ) == 0
+   !!
+   !!           ----------- WARNING -----------
+   !!           -------------------------------
+   !!
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE iom            !
+   USE ioipsl  , ONLY :   ymds2ju   ! for calendar
+   USE prtctl         ! Print control
+   USE trc_oce , ONLY : lk_offline ! offline flag
+   USE timing         ! Timing
+   USE restart        ! restart
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   day        ! called by step.F90
+   PUBLIC   day_init   ! called by istate.F90
+   PUBLIC   day_mth    ! Needed by TAM
+
+   INTEGER, PUBLIC ::   nsecd, nsecd05, ndt, ndt05   !: (PUBLIC for TAM)
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: daymod.F90 6140 2015-12-21 11:35:23Z timgraham $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE day_init
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE day_init  ***
+      !!
+      !! ** Purpose :   Initialization of the calendar values to their values 1 time step before nit000
+      !!                because day will be called at the beginning of step
+      !!
+      !! ** Action  : - nyear        : current year
+      !!              - nmonth       : current month of the year nyear
+      !!              - nday         : current day of the month nmonth
+      !!              - nday_year    : current day of the year nyear
+      !!              - nsec_year    : current time step counted in second since 00h jan 1st of the current year
+      !!              - nsec_month   : current time step counted in second since 00h 1st day of the current month
+      !!              - nsec_day     : current time step counted in second since 00h of the current day
+      !!              - nsec1jan000  : second since Jan. 1st 00h of nit000 year and Jan. 1st 00h of the current year
+      !!              - nmonth_len, nyear_len, nmonth_half, nmonth_end through day_mth
+      !!----------------------------------------------------------------------
+      INTEGER  ::   inbday, idweek
+      REAL(wp) ::   zjul
+      !!----------------------------------------------------------------------
+      !
+      ! max number of seconds between each restart
+      IF( REAL( nitend - nit000 + 1 ) * rdt > REAL( HUGE( nsec1jan000 ) ) ) THEN
+         CALL ctl_stop( 'The number of seconds between each restart exceeds the integer 4 max value: 2^31-1. ',   &
+            &           'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' )
+      ENDIF
+      ! all calendar staff is based on the fact that MOD( rday, rdt ) == 0
+      !SF IF( MOD( rday     , rdt   ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' )
+      !SF IF( MOD( rday     , 2.    ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
+      !SF IF( MOD( rdt      , 2.    ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
+      nsecd   = NINT(rday       )
+      nsecd05 = NINT(0.5 * rday )
+      ndt     = NINT(      rdt  )
+      ndt05   = NINT(0.5 * rdt  )
+
+      IF( .NOT. lk_offline ) CALL day_rst( nit000, 'READ' )
+
+      ! set the calandar from ndastp (read in restart file and namelist)
+
+      nyear   =   ndastp / 10000
+      nmonth  = ( ndastp - (nyear * 10000) ) / 100
+      nday    =   ndastp - (nyear * 10000) - ( nmonth * 100 )
+
+      nhour   =   nn_time0 / 100
+      nminute = ( nn_time0 - nhour * 100 )
+
+      CALL ymds2ju( nyear, nmonth, nday, nhour*3600._wp+nminute*60._wp, fjulday )  
+      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < 0.1 / rday )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
+      IF( nn_time0*3600 - ndt05 .lt. 0 ) fjulday = fjulday + 1.                    ! move back to the day at nit000 (and not at nit000 - 1)
+
+      nsec1jan000 = 0
+      CALL day_mth
+
+      IF ( nday == 0 ) THEN     !   for ex if ndastp = ndate0 - 1
+         nmonth = nmonth - 1
+         nday = nmonth_len(nmonth)
+      ENDIF
+      IF ( nmonth == 0 ) THEN   ! go at the end of previous year
+         nmonth = 12
+         nyear = nyear - 1
+         nsec1jan000 = nsec1jan000 - nsecd * nyear_len(0)
+         IF( nleapy == 1 )   CALL day_mth
+      ENDIF
+
+      ! day since january 1st
+      nday_year = nday + SUM( nmonth_len(1:nmonth - 1) )
+
+      !compute number of days between last monday and today
+      CALL ymds2ju( 1900, 01, 01, 0.0, zjul )  ! compute julian day value of 01.01.1900 (our reference that was a Monday)
+      inbday = FLOOR(fjulday - zjul)            ! compute nb day between  01.01.1900 and start of current day
+      idweek = MOD(inbday, 7)                  ! compute nb day between last monday and current day
+      IF (idweek .lt. 0) idweek=idweek+7       ! Avoid negative values for dates before 01.01.1900
+
+      ! number of seconds since the beginning of current year/month/week/day at the middle of the time-step
+      IF (nhour*3600+nminute*60-ndt05 .gt. 0) THEN
+         ! 1 timestep before current middle of first time step is still the same day
+         nsec_year  = (nday_year-1) * nsecd + nhour*3600+nminute*60 - ndt05 
+         nsec_month = (nday-1)      * nsecd + nhour*3600+nminute*60 - ndt05    
+      ELSE
+         ! 1 time step before the middle of the first time step is the previous day 
+         nsec_year  = nday_year * nsecd + nhour*3600+nminute*60 - ndt05 
+         nsec_month = nday      * nsecd + nhour*3600+nminute*60 - ndt05   
+      ENDIF
+      nsec_week  = idweek    * nsecd + nhour*3600+nminute*60 - ndt05
+      nsec_day   =             nhour*3600+nminute*60 - ndt05 
+      IF( nsec_day .lt. 0 ) nsec_day = nsec_day + nsecd
+      IF( nsec_week .lt. 0 ) nsec_week = nsec_week + nsecd*7
+
+      ! control print
+      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8,a,i8,a,i8)')' =======>> 1/2 time step before the start of the run DATE Y/M/D = ',   &
+           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week, '  &
+           &                   nsec_month:', nsec_month , '  nsec_year:' , nsec_year
+
+      ! Up to now, calendar parameters are related to the end of previous run (nit000-1)
+      ! call day to set the calendar parameters at the begining of the current simulaton. needed by iom_init
+      CALL day( nit000 )
+      !
+   END SUBROUTINE day_init
+
+
+   SUBROUTINE day_mth
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE day_init  ***
+      !!
+      !! ** Purpose :   calendar values related to the months
+      !!
+      !! ** Action  : - nmonth_len    : length in days of the months of the current year
+      !!              - nyear_len     : length in days of the previous/current year
+      !!              - nmonth_half   : second since the beginning of the year and the halft of the months
+      !!              - nmonth_end    : second since the beginning of the year and the end of the months
+      !!----------------------------------------------------------------------
+      INTEGER  ::   jm               ! dummy loop indice
+      !!----------------------------------------------------------------------
+
+      ! length of the month of the current year (from nleapy, read in namelist)
+      IF ( nleapy < 2 ) THEN
+         nmonth_len(:) = (/ 31, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31 /)
+         nyear_len(:) = 365
+         IF ( nleapy == 1 ) THEN   ! we are using calandar with leap years
+            IF ( MOD(nyear-1, 4) == 0 .AND. ( MOD(nyear-1, 400) == 0 .OR. MOD(nyear-1, 100) /= 0 ) ) THEN
+               nyear_len(0)  = 366
+            ENDIF
+            IF ( MOD(nyear  , 4) == 0 .AND. ( MOD(nyear  , 400) == 0 .OR. MOD(nyear  , 100) /= 0 ) ) THEN
+               nmonth_len(2) = 29
+               nyear_len(1)  = 366
+            ENDIF
+            IF ( MOD(nyear+1, 4) == 0 .AND. ( MOD(nyear+1, 400) == 0 .OR. MOD(nyear+1, 100) /= 0 ) ) THEN
+               nyear_len(2)  = 366
+            ENDIF
+         ENDIF
+      ELSE
+         nmonth_len(:) = nleapy   ! all months with nleapy days per year
+         nyear_len(:) = 12 * nleapy
+      ENDIF
+
+      ! half month in second since the begining of the year:
+      ! time since Jan 1st   0     1     2    ...    11    12    13
+      !          ---------*--|--*--|--*--| ... |--*--|--*--|--*--|--------------------------------------
+      !                 <---> <---> <--->  ...  <---> <---> <--->
+      ! month number      0     1     2    ...    11    12    13
+      !
+      ! nmonth_half(jm) = rday * REAL( 0.5 * nmonth_len(jm) + SUM(nmonth_len(1:jm-1)) )
+      nmonth_half(0) = - nsecd05 * nmonth_len(0)
+      DO jm = 1, 13
+         nmonth_half(jm) = nmonth_half(jm-1) + nsecd05 * ( nmonth_len(jm-1) + nmonth_len(jm) )
+      END DO
+
+      nmonth_end(0) = 0
+      DO jm = 1, 13
+         nmonth_end(jm) = nmonth_end(jm-1) + nsecd * nmonth_len(jm)
+      END DO
+      !
+   END SUBROUTINE
+
+
+   SUBROUTINE day( kt )
+      !!----------------------------------------------------------------------
+      !!                      ***  ROUTINE day  ***
+      !!
+      !! ** Purpose :   Compute the date with a day iteration IF necessary.
+      !!
+      !! ** Method  : - ???
+      !!
+      !! ** Action  : - nyear     : current year
+      !!              - nmonth    : current month of the year nyear
+      !!              - nday      : current day of the month nmonth
+      !!              - nday_year : current day of the year nyear
+      !!              - ndastp    : = nyear*10000 + nmonth*100 + nday
+      !!              - adatrj    : date in days since the beginning of the run
+      !!              - nsec_year : current time of the year (in second since 00h, jan 1st)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt        ! ocean time-step indices
+      !
+      CHARACTER (len=25) ::   charout
+      REAL(wp)           ::   zprec      ! fraction of day corresponding to 0.1 second
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('day')
+      !
+      zprec = 0.1 / rday
+      !                                                 ! New time-step
+      nsec_year  = nsec_year  + ndt
+      nsec_month = nsec_month + ndt
+      nsec_week  = nsec_week  + ndt
+      nsec_day   = nsec_day   + ndt
+      adatrj  = adatrj  + rdt / rday
+      fjulday = fjulday + rdt / rday
+      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < zprec )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
+      IF( ABS(adatrj  - REAL(NINT(adatrj ),wp)) < zprec )   adatrj  = REAL(NINT(adatrj ),wp)   ! avoid truncation error
+
+      IF( nsec_day > nsecd ) THEN                       ! New day
+         !
+         nday      = nday + 1
+         nday_year = nday_year + 1
+         nsec_day  = ndt05
+         !
+         IF( nday == nmonth_len(nmonth) + 1 ) THEN      ! New month
+            nday   = 1
+            nmonth = nmonth + 1
+            nsec_month = ndt05
+            IF( nmonth == 13 ) THEN                     ! New year
+               nyear     = nyear + 1
+               nmonth    = 1
+               nday_year = 1
+               nsec_year = ndt05
+               nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1)
+               IF( nleapy == 1 )   CALL day_mth
+            ENDIF
+         ENDIF
+         !
+         ndastp = nyear * 10000 + nmonth * 100 + nday   ! New date
+         !
+         !compute first day of the year in julian days
+         CALL ymds2ju( nyear, 01, 01, 0.0, fjulstartyear )
+         !
+         IF(lwp) WRITE(numout,'(a,i8,a,i4.4,a,i2.2,a,i2.2,a,i3.3)') '======>> time-step =', kt,   &
+              &   '      New day, DATE Y/M/D = ', nyear, '/', nmonth, '/', nday, '      nday_year = ', nday_year
+         IF(lwp) WRITE(numout,'(a,i8,a,i7,a,i5)') '         nsec_year = ', nsec_year,   &
+              &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day, '   nsec_week = ', nsec_week
+      ENDIF
+
+      IF( nsec_week > 7*nsecd )   nsec_week = ndt05     ! New week
+
+      IF(ln_ctl) THEN
+         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
+         CALL prt_ctl_info(charout)
+      ENDIF
+
+      IF( .NOT. lk_offline ) CALL rst_opn( kt )               ! Open the restart file if needed and control lrst_oce
+      IF( lrst_oce         ) CALL day_rst( kt, 'WRITE' )      ! write day restart information
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('day')
+      !
+   END SUBROUTINE day
+
+
+   SUBROUTINE day_rst( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ts_rst  ***
+      !!
+      !!  ** Purpose : Read or write calendar in restart file:
+      !!
+      !!  WRITE(READ) mode:
+      !!       kt        : number of time step since the begining of the experiment at the
+      !!                   end of the current(previous) run
+      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
+      !!                   end of the current(previous) run (REAL -> keep fractions of day)
+      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
+      !!
+      !!   According to namelist parameter nrstdt,
+      !!       nrstdt = 0  no control on the date (nit000 is  arbitrary).
+      !!       nrstdt = 1  we verify that nit000 is equal to the last
+      !!                   time step of previous run + 1.
+      !!       In both those options, the  exact duration of the experiment
+      !!       since the beginning (cumulated duration of all previous restart runs)
+      !!       is not stored in the restart and is assumed to be (nit000-1)*rdt.
+      !!       This is valid is the time step has remained constant.
+      !!
+      !!       nrstdt = 2  the duration of the experiment in days (adatrj)
+      !!                    has been stored in the restart file.
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      !
+      REAL(wp) ::   zkt, zndastp, zdayfrac, ksecs, ktime
+      INTEGER  ::   ihour, iminute
+      !!----------------------------------------------------------------------
+
+      IF( TRIM(cdrw) == 'READ' ) THEN
+
+         IF( iom_varid( numror, 'kt', ldstop = .FALSE. ) > 0 ) THEN
+            ! Get Calendar informations
+            CALL iom_get( numror, 'kt', zkt )   ! last time-step of previous run
+            IF(lwp) THEN
+               WRITE(numout,*) ' *** Info read in restart : '
+               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
+               WRITE(numout,*) ' *** restart option'
+               SELECT CASE ( nrstdt )
+               CASE ( 0 )   ;   WRITE(numout,*) ' nrstdt = 0 : no control of nit000'
+               CASE ( 1 )   ;   WRITE(numout,*) ' nrstdt = 1 : no control the date at nit000 (use ndate0 read in the namelist)'
+               CASE ( 2 )   ;   WRITE(numout,*) ' nrstdt = 2 : calendar parameters read in restart'
+               END SELECT
+               WRITE(numout,*)
+            ENDIF
+            ! Control of date
+            IF( nit000 - NINT( zkt ) /= 1 .AND. nrstdt /= 0 )                                         &
+                 &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart',                 &
+                 &                  ' verify the restart file or rerun with nrstdt = 0 (namelist)' )
+            ! define ndastp and adatrj
+            IF ( nrstdt == 2 ) THEN
+               ! read the parameters corresponding to nit000 - 1 (last time step of previous run)
+               CALL iom_get( numror, 'ndastp', zndastp )
+               ndastp = NINT( zndastp )
+               CALL iom_get( numror, 'adatrj', adatrj  )
+	       CALL iom_get( numror, 'ntime', ktime )
+	       nn_time0=INT(ktime)
+               ! calculate start time in hours and minutes
+	       zdayfrac=adatrj-INT(adatrj)
+	       ksecs = NINT(zdayfrac*86400)	       ! Nearest second to catch rounding errors in adatrj	       
+	       ihour = INT(ksecs/3600)
+	       iminute = ksecs/60-ihour*60
+ 	        
+               ! Add to nn_time0
+               nhour   =   nn_time0 / 100
+               nminute = ( nn_time0 - nhour * 100 )
+	       nminute=nminute+iminute
+	       
+	       IF( nminute >= 60 ) THEN
+	          nminute=nminute-60
+		  nhour=nhour+1
+	       ENDIF
+	       nhour=nhour+ihour
+	       IF( nhour >= 24 ) THEN
+		  nhour=nhour-24
+	          adatrj=adatrj+1
+	       ENDIF	       
+	       nn_time0 = nhour * 100 + nminute
+	       adatrj = INT(adatrj)                    ! adatrj set to integer as nn_time0 updated	       
+            ELSE
+               ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day)
+               ndastp = ndate0        ! ndate0 read in the namelist in dom_nam
+               nhour   =   nn_time0 / 100
+               nminute = ( nn_time0 - nhour * 100 )
+               IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0)
+               adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday
+               ! note this is wrong if time step has changed during run
+            ENDIF
+         ELSE
+            ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day)
+            ndastp = ndate0           ! ndate0 read in the namelist in dom_nam
+            nhour   =   nn_time0 / 100
+	    nminute = ( nn_time0 - nhour * 100 )
+            IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0)
+            adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday
+         ENDIF
+         IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) ' *** Info used values : '
+            WRITE(numout,*) '   date ndastp                                      : ', ndastp
+            WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
+	    WRITE(numout,*) '   nn_time0                                         : ',nn_time0
+            WRITE(numout,*)
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+         !
+         IF( kt == nitrst ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'rst_write : write oce restart file  kt =', kt
+            IF(lwp) WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         ! calendar control
+         CALL iom_rstput( kt, nitrst, numrow, 'kt'     , REAL( kt    , wp) )   ! time-step
+         CALL iom_rstput( kt, nitrst, numrow, 'ndastp' , REAL( ndastp, wp) )   ! date
+         CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj            )   ! number of elapsed days since
+         !                                                                     ! the begining of the run [s]
+	 CALL iom_rstput( kt, nitrst, numrow, 'ntime'  , REAL( nn_time0, wp) ) ! time
+      ENDIF
+      !
+   END SUBROUTINE day_rst
+
+   !!======================================================================
+END MODULE daymod
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/traadv_fct.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/traadv_fct.F90	(revision 6879)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/traadv_fct.F90	(revision 6879)
@@ -0,0 +1,795 @@
+MODULE traadv_fct
+   !!==============================================================================
+   !!                       ***  MODULE  traadv_fct  ***
+   !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
+   !!==============================================================================
+   !! History :  3.7  !  2015-09  (L. Debreu, G. Madec)  original code (inspired from traadv_tvd.F90)
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  tra_adv_fct    : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
+   !!  tra_adv_fct_zts: update the tracer trend with a 3D advective trends using a 2nd order FCT scheme 
+   !!                   with sub-time-stepping in the vertical direction
+   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm 
+   !!  interp_4th_cpt : 4th order compact scheme for the vertical component of the advection
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE trc_oce        ! share passive tracers/Ocean variables
+   USE trd_oce        ! trends: ocean variables
+   USE trdtra         ! tracers trends
+   USE diaptr         ! poleward transport diagnostics
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_adv_fct        ! routine called by traadv.F90
+   PUBLIC   tra_adv_fct_zts    ! routine called by traadv.F90
+   PUBLIC   interp_4th_cpt     ! routine called by traadv_cen.F90
+
+   LOGICAL  ::   l_trd   ! flag to compute trends
+   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! $Id: traadv_fct.F90 6140 2015-12-21 11:35:23Z timgraham $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn,       &
+      &                                              ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_fct  ***
+      !! 
+      !! **  Purpose :   Compute the now trend due to total advection of tracers
+      !!               and add it to the general trend of tracer equations
+      !!
+      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
+      !!               (choice through the value of kn_fct)
+      !!               - on the vertical the 4th order is a compact scheme 
+      !!               - corrected flux (monotonic correction) 
+      !!
+      !! ** Action : - update pta  with the now advective tracer trends
+      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
+      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
+      !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
+      INTEGER                              , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
+      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
+      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
+      !
+      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices  
+      REAL(wp) ::   ztra                                     ! local scalar
+      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      -
+      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      l_trd = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
+      !
+      IF( l_trd )  THEN
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
+         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
+      ENDIF
+      !
+      !                          ! surface & bottom value : flux set to zero one for all
+      zwz(:,:, 1 ) = 0._wp            
+      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp
+      !
+      zwi(:,:,:) = 0._wp        
+      !
+      DO jn = 1, kjpt            !==  loop over the tracers  ==!
+         !
+         !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
+         !                    !* upstream tracer flux in the i and j direction 
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ! upstream scheme
+                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
+                  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) )
+                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
+               END DO
+            END DO
+         END DO
+         !                    !* upstream tracer flux in the k direction *!
+         DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+            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) )
+                  zwz(ji,jj,jk) = 0.5 * ( 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( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz 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
+                     zwz(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
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+            ENDIF
+         ENDIF
+         !               
+         DO jk = 1, jpkm1     !* trend and after field with monotonic scheme
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! total intermediate advective trends
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  ! update and guess with monotonic sheme
+!!gm why tmask added in the two following lines ???    the mask is done in tranxt !
+                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)        + ztra   * tmask(ji,jj,jk)
+!!gm BUG                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk)
+                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + e3t_n(ji,jj,jk) * p2dt * ztra )   &
+                     &             /   e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( zwi, 'T', 1. )  ! Lateral boundary conditions on zwi  (unchanged sign)
+         !                
+         IF( l_trd )  THEN             ! trend diagnostics (contribution of upstream fluxes)
+            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
+         END IF
+         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
+         ENDIF
+         !
+         !        !==  anti-diffusive flux : high order minus low order  ==!
+         !
+         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes
+         !
+         CASE(  2  )                   !- 2nd order centered
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+         CASE(  4  )                   !- 4th order centered
+            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero
+            zltv(:,:,jpk) = 0._wp
+            DO jk = 1, jpkm1                 ! Laplacian
+               DO jj = 1, jpjm1                    ! 1st derivative (gradient)
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                  END DO
+               END DO
+               DO jj = 2, jpjm1                    ! 2nd derivative * 1/ 6
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
+                     zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO jk = 1, jpkm1                 ! Horizontal advective fluxes
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points
+                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
+                     !                                                  ! C4 minus upstream advective fluxes 
+                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
+                  END DO
+               END DO
+            END DO         
+            !
+         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
+            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
+            ztv(:,:,jpk) = 0._wp
+            DO jk = 1, jpkm1                 ! 1st derivative (gradient)
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
+            !
+            DO jk = 1, jpkm1                 ! Horizontal advective fluxes
+               DO jj = 2, jpjm1
+                  DO ji = 2, fs_jpim1   ! vector opt.
+                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points (x2)
+                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
+                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
+                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) )
+                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) )
+                     !                                                  ! C4 minus upstream advective fluxes 
+                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+         END SELECT
+         !                      
+         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values)
+         !
+         CASE(  2  )                   !- 2nd order centered
+            DO jk = 2, jpkm1    
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
+                        &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+         CASE(  4  )                   !- 4th order COMPACT
+            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point
+            DO jk = 2, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+         END SELECT
+         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
+            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
+         ENDIF
+         !
+         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
+         CALL lbc_lnk( zwz, 'W',  1. )
+         !
+         !        !==  monotonicity algorithm  ==!
+         !
+         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
+         !
+         !        !==  final trend with corrected fluxes  ==!
+         !
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.  
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) &
+                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF( l_trd ) THEN     ! trend diagnostics (contribution of upstream fluxes)
+            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
+            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
+            !
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
+            !
+            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
+         END IF
+         !                    ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )   htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) )
+           IF( jn == jp_sal )   str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) )
+         ENDIF
+         !
+      END DO                     ! end of tracer loop
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct')
+      !
+   END SUBROUTINE tra_adv_fct
+
+
+   SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
+      &                                                  ptb, ptn, pta, kjpt, kn_fct_zts )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_fct_zts  ***
+      !! 
+      !! **  Purpose :   Compute the now trend due to total advection of 
+      !!       tracers and add it to the general trend of tracer equations
+      !!
+      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with
+      !!       corrected flux (monotonic correction). This version use sub-
+      !!       timestepping for the vertical advection which increases stability
+      !!       when vertical metrics are small.
+      !!       note: - this advection scheme needs a leap-frog time scheme
+      !!
+      !! ** Action : - update (pta) with the now advective tracer trends
+      !!             - save the trends 
+      !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
+      INTEGER                              , INTENT(in   ) ::   kn_fct_zts      ! number of number of vertical sub-timesteps
+      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
+      !
+      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection
+      REAL(wp)                                             ::   zr_p2dt         ! reciprocal of tracer timestep
+      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices  
+      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
+      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps
+      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection
+      REAL(wp) ::   ztra            ! local scalar
+      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
+      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zwx_sav , zwy_sav
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdx, ztrdy, ztrdz
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct_zts')
+      !
+      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
+      CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
+      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      l_trd = .FALSE.
+      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
+      !
+      IF( l_trd )  THEN
+         CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
+         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
+      ENDIF
+      !
+      zwi(:,:,:) = 0._wp
+      z_rzts = 1._wp / REAL( kn_fct_zts, wp )
+      zr_p2dt = 1._wp / p2dt
+      !
+      ! surface & Bottom value : flux set to zero for all tracers
+      zwz(:,:, 1 ) = 0._wp
+      zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp
+      zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp
+      !
+      !                                                          ! ===========
+      DO jn = 1, kjpt                                            ! tracer loop
+         !                                                       ! ===========
+         !
+         ! Upstream advection with initial mass fluxes & intermediate update
+         DO jk = 1, jpkm1        ! upstream tracer flux in the i and j direction
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ! upstream scheme
+                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
+                  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) )
+                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
+               END DO
+            END DO
+         END DO
+         !                       ! upstream tracer flux in the k direction
+         DO jk = 2, jpkm1              ! Interior value
+            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) )
+                  zwz(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( ln_linssh ) THEN          ! top value : linear free surface case only (as zwz is multiplied by wmask)
+            IF( ln_isfcav ) THEN             ! ice-shelf cavities: top value
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 
+                  END DO
+               END DO   
+            ELSE                             ! no cavities, surface value
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+            ENDIF
+         ENDIF
+         !
+         DO jk = 1, jpkm1         ! total advective trend
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !                             ! total intermediate advective trends
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  !                             ! update and guess with monotonic sheme
+                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra
+                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !                           
+         CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign)
+         !                
+         IF( l_trd )  THEN                ! trend diagnostics (contribution of upstream fluxes)
+            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
+         END IF
+         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
+         ENDIF
+
+         ! 3. anti-diffusive flux : high order minus low order
+         ! ---------------------------------------------------
+
+         DO jk = 1, jpkm1                    !* horizontal anti-diffusive fluxes
+            !
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwx_sav(ji,jj) = zwx(ji,jj,jk)
+                  zwy_sav(ji,jj) = zwy(ji,jj,jk)
+                  !
+                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
+                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
+               END DO
+            END DO
+            !
+            DO jj = 2, jpjm1                    ! partial horizontal divergence
+               DO ji = fs_2, fs_jpim1
+                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
+                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
+               END DO
+            END DO
+            !
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj)
+                  zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj)
+               END DO
+            END DO
+         END DO
+         !
+         !                                !* vertical anti-diffusive flux
+         zwz_sav(:,:,:)   = zwz(:,:,:)
+         ztrs   (:,:,:,1) = ptb(:,:,:,jn)
+         zwzts  (:,:,:)   = 0._wp
+         !
+         DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop
+            !
+            IF( jl == 1 ) THEN                        ! Euler forward to kick things off
+               jtb = 1   ;   jtn = 1   ;   jta = 2
+               zts(:) = p2dt * z_rzts
+               jtaken = MOD( kn_fct_zts + 1 , 2)            ! Toggle to collect every second flux
+               !                                            ! starting at jl =1 if kn_fct_zts is odd; 
+               !                                            ! starting at jl =2 otherwise
+            ELSEIF( jl == 2 ) THEN                    ! First leapfrog step
+               jtb = 1   ;   jtn = 2   ;   jta = 3
+               zts(:) = 2._wp * p2dt * z_rzts
+            ELSE                                      ! Shuffle pointers for subsequent leapfrog steps
+               jtb = MOD(jtb,3) + 1
+               jtn = MOD(jtn,3) + 1
+               jta = MOD(jta,3) + 1
+            ENDIF
+            DO jk = 2, jpkm1                          ! interior value
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) * wmask(ji,jj,jk)
+                     IF( jtaken == 0 )   zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk)    ! Accumulate time-weighted vertcal flux
+                  END DO
+               END DO
+            END DO
+            IF( ln_linssh ) THEN                    ! top value (only in linear free surface case)
+               IF( ln_isfcav ) THEN                      ! ice-shelf cavities
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        zwz(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 ocean cavities
+                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+               ENDIF
+            ENDIF
+            !
+            jtaken = MOD( jtaken + 1 , 2 )
+            !
+            DO jk = 2, jpkm1                             ! total advective trends
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb)                                                 &
+                        &               - zts(jk) * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
+                        &                         * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+         END DO
+
+         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
+         CALL lbc_lnk( zwz, 'W',  1. )
+
+         ! 4. monotonicity algorithm
+         ! -------------------------
+         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
+
+
+         ! 5. final trend with corrected fluxes
+         ! ------------------------------------
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.  
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )       &
+                     &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   )   &
+                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+
+         !                                 ! trend diagnostics (contribution of upstream fluxes)
+         IF( l_trd )  THEN 
+            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
+            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
+            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
+            !
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )   
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )  
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
+            !
+            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
+         END IF
+         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:)
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:)
+         ENDIF
+         !
+      END DO
+      !
+      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
+      CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
+      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts')
+      !
+   END SUBROUTINE tra_adv_fct_zts
+
+
+   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE nonosc  ***
+      !!     
+      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
+      !!       scheme and the before field by a nonoscillatory algorithm 
+      !!
+      !! **  Method  :   ... ???
+      !!       warning : pbef and paft must be masked, but the boundaries
+      !!       conditions on the fluxes are not necessary zalezak (1979)
+      !!       drange (1995) multi-dimensional forward-in-time and upstream-
+      !!       in-space based differencing for fluid
+      !!----------------------------------------------------------------------
+      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
+      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ikm1         ! local integer
+      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
+      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('nonosc')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
+      !
+      zbig  = 1.e+40_wp
+      zrtrn = 1.e-15_wp
+      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
+
+      ! Search local extrema
+      ! --------------------
+      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
+      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   &
+         &        paft * tmask - zbig * ( 1._wp - tmask )  )
+      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   &
+         &        paft * tmask + zbig * ( 1._wp - tmask )  )
+
+      DO jk = 1, jpkm1
+         ikm1 = MAX(jk-1,1)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+
+               ! search maximum in neighbourhood
+               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
+                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
+                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
+                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
+
+               ! search minimum in neighbourhood
+               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
+                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
+                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
+                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
+
+               ! positive part of the flux
+               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
+                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
+                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
+
+               ! negative part of the flux
+               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
+                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
+                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
+
+               ! up & down beta terms
+               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
+               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
+               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
+
+      ! 3. monotonic flux in the i & j direction (paa & pbb)
+      ! ----------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
+               zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
+               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
+               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
+
+               zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
+               zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
+               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
+               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
+
+      ! monotonic flux in the k direction, i.e. pcc
+      ! -------------------------------------------
+               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
+               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
+               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
+               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('nonosc')
+      !
+   END SUBROUTINE nonosc
+
+
+   SUBROUTINE interp_4th_cpt( pt_in, pt_out )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE interp_4th_cpt  ***
+      !! 
+      !! **  Purpose :   Compute the interpolation of tracer at w-point
+      !!
+      !! **  Method  :   4th order compact interpolation
+      !!----------------------------------------------------------------------
+      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! now tracer fields
+      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! now tracer field interpolated at w-pts
+      !
+      INTEGER :: ji, jj, jk   ! dummy loop integers
+      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
+      !!----------------------------------------------------------------------
+      
+      DO jk = 3, jpkm1        !==  build the three diagonal matrix  ==!
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zwd (ji,jj,jk) = 4._wp
+               zwi (ji,jj,jk) = 1._wp
+               zws (ji,jj,jk) = 1._wp
+               zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
+               !
+               IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom
+                  zwd (ji,jj,jk) = 1._wp
+                  zwi (ji,jj,jk) = 0._wp
+                  zws (ji,jj,jk) = 0._wp
+                  zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )    
+               ENDIF
+            END DO
+         END DO
+      END DO
+      !
+      jk=2                                            ! Switch to second order centered at top
+      DO jj=1,jpj
+         DO ji=1,jpi
+            zwd (ji,jj,jk) = 1._wp
+            zwi (ji,jj,jk) = 0._wp
+            zws (ji,jj,jk) = 0._wp
+            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
+         END DO
+      END DO   
+      !
+      !                       !==  tridiagonal solve  ==!
+      DO jj = 1, jpj                ! first recurrence
+         DO ji = 1, jpi
+            zwt(ji,jj,2) = zwd(ji,jj,2)
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 1, jpj                ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+         DO ji = 1, jpi
+            pt_out(ji,jj,2) = zwrm(ji,jj,2)
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
+            END DO
+         END DO
+      END DO
+
+      DO jj = 1, jpj                ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
+         DO ji = 1, jpi
+            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 2, -1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !    
+   END SUBROUTINE interp_4th_cpt
+   
+   !!======================================================================
+END MODULE traadv_fct
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/usrdef_sbc.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/usrdef_sbc.F90	(revision 6879)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/usrdef_sbc.F90	(revision 6879)
@@ -0,0 +1,145 @@
+MODULE usrdef_sbc
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_sbc  ***
+   !! Ocean forcing:  user defined momentum, heat and freshwater forcings
+   !!
+   !!                ===     Here  GYRE configuration      ===
+   !!
+   !!=====================================================================
+   !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)  user defined interface
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_sbc    : user defined surface bounday conditions in GYRE case
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE phycst          ! physical constants
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_fortran     !
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_sbc    ! routine called in sbcmod module
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_sbc( kt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE usr_def_sbc  ***
+      !!              
+      !! ** Purpose :   provide at each time-step the surface boundary
+      !!              condition, i.e. the momentum, heat and freshwater fluxes.
+      !!
+      !! ** Method  :   analytical seasonal cycle (all 0 fields, for overflow
+      !!              case).
+      !!                CAUTION : never mask the surface stress field !
+      !!
+      !! ** Action  : - set the ocean surface boundary condition, i.e.   
+      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!
+      !SF INTEGER  ::   ji, jj                 ! dummy loop indices
+      !SF INTEGER  ::   zyear0                 ! initial year 
+      !SF INTEGER  ::   zmonth0                ! initial month
+      !SF INTEGER  ::   zday0                  ! initial day
+      !SF INTEGER  ::   zday_year0             ! initial day since january 1st
+      !SF REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
+      !SF REAL(wp) ::   ztime                  ! time in hour
+      !SF REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
+      !SF REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
+      !SF REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
+      !SF REAL(wp) ::   ztaun                  ! intensity
+      !SF REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
+      !SF REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
+      !SF REAL(wp) ::   zsumemp, zsurf
+      !SF REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
+      !SF REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
+      !SF REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
+      !SF REAL(wp) ::   zyydd                 ! number of days in one year
+!SF -------
+!SF sbc ana:
+      REAL(wp) ::   zrhoa  = 1.22_wp      ! air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3_wp    ! drag coefficient
+      REAL(wp) ::   zfact, ztx            ! local scalars
+      REAL(wp) ::   zcoef, zty, zmod      !   -      -
+      INTEGER  ::   nn_tau000   =   0     !  gently increase the stress over the first ntau_rst time-steps
+      REAL(wp) ::   rn_utau0    =   0.e0  !  uniform value for the i-stress
+      REAL(wp) ::   rn_vtau0    =   0.e0  !  uniform value for the j-stress
+      REAL(wp) ::   rn_qns0     =   0.e0  !  uniform value for the total heat flux
+      REAL(wp) ::   rn_qsr0     =   0.e0  !  uniform value for the solar radiation
+      REAL(wp) ::   rn_emp0     =   0.e0  !  uniform value for the freswater budget (E-P) 
+      !
+      !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF(lwp) WRITE(numout,*)' sbc_ana : Constant surface fluxes read in namsbc_ana namelist'
+         IF(lwp) WRITE(numout,*)' ~~~~~~~ '
+         IF(lwp) WRITE(numout,*)'              spin up of the stress  nn_tau000 = ', nn_tau000, ' time-steps'
+         IF(lwp) WRITE(numout,*)'              constant j-stress      rn_vtau0  = ', rn_vtau0 , ' N/m2'
+         IF(lwp) WRITE(numout,*)'              non solar heat flux    rn_qns0 = ', rn_qns0  , ' W/m2'
+         IF(lwp) WRITE(numout,*)'              solar heat flux        rn_qsr0 = ', rn_qsr0  , ' W/m2'
+         IF(lwp) WRITE(numout,*)'              net heat flux          rn_emp0 = ', rn_emp0  , ' Kg/m2/s'
+         !
+         nn_tau000 = MAX( nn_tau000, 1 )     ! must be >= 1
+         !
+         utau(:,:) = rn_utau0
+         vtau(:,:) = rn_vtau0
+         taum(:,:) = SQRT ( rn_utau0 * rn_utau0 + rn_vtau0 * rn_vtau0 )
+         wndm(:,:) = SQRT ( taum(1,1) /  ( zrhoa * zcdrag ) )
+         !
+         emp (:,:) = rn_emp0
+         sfx (:,:) = 0.0_wp
+         qns (:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp      ! including heat content associated with mass flux at SST
+         qsr (:,:) = rn_qsr0
+         !         
+      ENDIF
+
+      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+         !
+         IF( kt <= nn_tau000 ) THEN       ! Increase the stress to its nominal value 
+            !                             ! during the first nn_tau000
+            !                             time-steps
+            zfact = 0.5 * (  1. - COS( rpi * REAL( kt, wp ) / REAL( nn_tau000, wp ) )  )
+            zcoef = 1. / ( zrhoa * zcdrag )
+            ztx   = zfact * rn_utau0
+            zty   = zfact * rn_vtau0
+            zmod  = SQRT( ztx * ztx + zty * zty )
+            utau(:,:) = ztx
+            vtau(:,:) = zty
+            taum(:,:) = zmod
+            zmod = SQRT( zmod * zcoef )
+            wndm(:,:) = zmod
+         ENDIF
+         !                                ! update heat and fresh water fluxes
+         !                                ! as they may have been changed by
+         !                                sbcssr module
+         emp (:,:) = rn_emp0              ! NB: qns changes with SST if emp /= 0
+         sfx (:,:) = 0._wp
+         qns (:,:) = rn_qns0 - emp(:,:) * sst_m(:,:) * rcp
+         qsr (:,:) = rn_qsr0
+         !
+      ENDIF
+      !
+!SF END sbc ana
+!SF -----------
+
+   END SUBROUTINE usr_def_sbc
+
+   !!======================================================================
+END MODULE usrdef_sbc
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/zdfini.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/zdfini.F90	(revision 6879)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/zdfini.F90	(revision 6879)
@@ -0,0 +1,179 @@
+MODULE zdfini
+   !!======================================================================
+   !!                      ***  MODULE  zdfini  ***
+   !! Ocean physics :   read vertical mixing namelist and check consistancy
+   !!======================================================================
+   !! History :  8.0  ! 1997-06  (G. Madec)  Original code from inimix
+   !!            1.0  ! 2002-08  (G. Madec)  F90 : free form
+   !!             -   ! 2005-06  (C. Ethe) KPP scheme
+   !!             -   ! 2009-07  (G. Madec) add avmb, avtb in restart for cen2 advection
+   !!            3.7  ! 2014-12  (G. Madec) remove KPP scheme
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   zdf_init    : initialization, namelist read, and parameters control
+   !!----------------------------------------------------------------------
+   USE par_oce         ! mesh and scale factors
+   USE zdf_oce         ! TKE vertical mixing          
+   USE sbc_oce         ! surface module (only for nn_isf in the option compatibility test)
+   USE zdftke          ! TKE vertical mixing
+   USE zdfgls          ! GLS vertical mixing
+   USE zdfric          ! Richardson vertical mixing   
+   USE zdfddm          ! double diffusion mixing      
+   USE zdfevd          ! enhanced vertical diffusion  
+   USE tranpc          ! convection: non penetrative adjustment
+   USE ldfslp          ! iso-neutral slopes
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! IOM library
+   USE lib_mpp         ! distribued memory computing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_init   ! routine called by opa.F90
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: zdfini.F90 5836 2015-10-26 14:49:40Z cetlod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE zdf_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_init  ***
+      !! 
+      !! ** Purpose :   initializations of the vertical ocean physics
+      !!
+      !! ** Method  :   Read namelist namzdf, control logicals 
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio, ios       ! local integers
+      !!
+      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,   &
+         &              ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namzdf in reference namelist : Vertical mixing parameters
+      READ  ( numnam_ref, namzdf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namzdf in reference namelist : Vertical mixing parameters
+      READ  ( numnam_cfg, namzdf, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzdf )
+
+      IF(lwp) THEN               !* Parameter print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_init: vertical physics'
+         WRITE(numout,*) '~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters'
+         WRITE(numout,*) '      vertical eddy viscosity             rn_avm0   = ', rn_avm0
+         WRITE(numout,*) '      vertical eddy diffusivity           rn_avt0   = ', rn_avt0
+         WRITE(numout,*) '      constant background or profile      nn_avb    = ', nn_avb
+         WRITE(numout,*) '      horizontal variation for avtb       nn_havtb  = ', nn_havtb
+         WRITE(numout,*) '      time splitting / backward scheme    ln_zdfexp = ', ln_zdfexp
+         WRITE(numout,*) '      number of time step                 nn_zdfexp = ', nn_zdfexp
+         WRITE(numout,*) '      enhanced vertical diffusion         ln_zdfevd = ', ln_zdfevd
+         WRITE(numout,*) '         applied on momentum (=1/0)       nn_evdm   = ', nn_evdm
+         WRITE(numout,*) '      vertical coefficient for evd        rn_avevd  = ', rn_avevd
+         WRITE(numout,*) '      non-penetrative convection (npc)    ln_zdfnpc = ', ln_zdfnpc
+         WRITE(numout,*) '      npc call  frequency                 nn_npc    = ', nn_npc
+         WRITE(numout,*) '      npc print frequency                 nn_npcp   = ', nn_npcp
+      ENDIF
+
+      !                          !* Parameter & logical controls
+      !                          !  ----------------------------
+      !
+      !                               ! ... check of vertical mixing scheme on tracers
+      !                                              ==> will be done in trazdf module
+      !
+      !                               ! ... check of mixing coefficient
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '   vertical mixing option :'
+      ioptio = 0
+      IF( lk_zdfcst ) THEN
+         IF(lwp) WRITE(numout,*) '      constant eddy diffusion coefficients'
+         ioptio = ioptio+1
+      ENDIF
+      IF( lk_zdfric ) THEN
+         IF(lwp) WRITE(numout,*) '      Richardson dependent eddy coefficients'
+         ioptio = ioptio+1
+      ENDIF
+      IF( lk_zdftke ) THEN
+         IF(lwp) WRITE(numout,*) '      TKE dependent eddy coefficients'
+         ioptio = ioptio+1
+      ENDIF
+      IF( lk_zdfgls ) THEN
+         IF(lwp) WRITE(numout,*) '      GLS dependent eddy coefficients'
+         ioptio = ioptio+1
+      ENDIF
+      IF( ioptio == 0 .OR. ioptio > 1 )   &
+         &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' )
+      IF( ( lk_zdfric .OR. lk_zdfgls ) .AND. ln_isfcav )   &
+         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' )
+      !
+      !                               ! ... Convection
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '   convection :'
+      !
+#if defined key_top
+      IF( ln_zdfnpc )   CALL ctl_stop( ' zdf_init: npc scheme is not working with key_top' )
+#endif
+      !
+      ioptio = 0
+      IF( ln_zdfnpc ) THEN
+         IF(lwp) WRITE(numout,*) '      use non penetrative convective scheme'
+         ioptio = ioptio+1
+      ENDIF
+      IF( ln_zdfevd ) THEN
+         IF(lwp) WRITE(numout,*) '      use enhanced vertical dif. scheme'
+         ioptio = ioptio+1
+      ENDIF
+      IF( lk_zdftke ) THEN
+         IF(lwp) WRITE(numout,*) '      use the 1.5 turbulent closure'
+      ENDIF
+      IF( lk_zdfgls ) THEN
+         IF(lwp) WRITE(numout,*) '      use the GLS closure scheme'
+      ENDIF
+      IF ( ioptio > 1 )   CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' )
+      IF( ioptio == 0 .AND. .NOT.( lk_zdftke .OR. lk_zdfgls ) )           &
+         !SF CALL ctl_stop( ' except for TKE or GLS physics, a convection scheme is',   &
+         CALL ctl_warn( ' except for TKE or GLS physics, a convection scheme is',   &
+         &              ' required: ln_zdfevd or ln_zdfnpc logicals' )
+
+      !                               !* Background eddy viscosity and diffusivity profil
+      IF( nn_avb == 0 ) THEN                ! Define avmb, avtb from namelist parameter
+         avmb(:) = rn_avm0
+         avtb(:) = rn_avt0                     
+      ELSE                                  ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
+         avmb(:) = rn_avm0
+         avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:)   ! m2/s
+         IF(ln_sco .AND. lwp)   CALL ctl_warn( 'avtb profile not valid in sco' )
+      ENDIF
+      !
+      IF( ln_rstart ) THEN                  !  Read avmb, avtb in restart (if exist)
+         ! if ln_traadv_cen, avmb, avtb have been modified in traadv_cen2 module. 
+         ! To ensure the restartability, avmb & avtb are written in the restart 
+         ! file in traadv_cen2 end read here. 
+         IF( iom_varid( numror, 'avmb', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numror, jpdom_unknown, 'avmb', avmb )
+            CALL iom_get( numror, jpdom_unknown, 'avtb', avtb )
+         ENDIF
+      ENDIF
+      !                                     ! 2D shape of the avtb
+      avtb_2d(:,:) = 1.e0                        ! uniform 
+      !
+      IF( nn_havtb == 1 ) THEN                   ! decrease avtb in the equatorial band
+           !  -15S -5S : linear decrease from avt0 to avt0/10.
+           !  -5S  +5N : cst value avt0/10.
+           !   5N  15N : linear increase from avt0/10, to avt0
+           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
+           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
+           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
+      ENDIF
+      !
+   END SUBROUTINE zdf_init
+
+   !!======================================================================
+END MODULE zdfini
