Changeset 5770

Timestamp:

2015-10-01T14:48:08+02:00 (9 years ago)

Author:

gm

Message:

#1593: LDF-ADV, step II.2: phasing the improvements/simplifications of advective tracer trend (see wiki page)

Location:

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO

Files:

• OPA_SRC/DOM/dom_oce.F90 (modified) (2 diffs)
• OPA_SRC/DOM/domain.F90 (modified) (4 diffs)
• OPA_SRC/DOM/dommsk.F90 (modified) (4 diffs)
• OPA_SRC/DOM/domngb.F90 (modified) (3 diffs)
• OPA_SRC/DOM/domwri.F90 (modified) (1 diff)
• OPA_SRC/DOM/domzgr.F90 (modified) (2 diffs)
• OPA_SRC/DOM/dtatsd.F90 (modified) (1 diff)
• OPA_SRC/DOM/istate.F90 (modified) (4 diffs)
• OPA_SRC/DYN/divcur.F90 (modified) (6 diffs)
• OPA_SRC/DYN/dynspg.F90 (modified) (1 diff)
• OPA_SRC/DYN/dynspg_flt.F90 (modified) (3 diffs)
• OPA_SRC/DYN/sshwzv.F90 (modified) (6 diffs)
• OPA_SRC/LBC/cla.F90 (deleted)
• OPA_SRC/TRA/traadv.F90 (modified) (8 diffs)
• OPA_SRC/TRA/traadv_cen.F90 (added)
• OPA_SRC/TRA/traadv_cen2.F90 (deleted)
• OPA_SRC/TRA/traadv_fct.F90 (moved) (moved from branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90) (21 diffs)
• OPA_SRC/TRA/traadv_mus.F90 (moved) (moved from branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90) (12 diffs)
• OPA_SRC/TRA/traadv_muscl2.F90 (deleted)
• OPA_SRC/TRA/traadv_qck.F90 (modified) (5 diffs)
• OPA_SRC/TRA/traadv_ubs.F90 (modified) (19 diffs)
• OPA_SRC/TRA/traldf.F90 (modified) (5 diffs)
• OPA_SRC/TRA/traldf_blp.F90 (modified) (5 diffs)
• OPA_SRC/TRA/trazdf.F90 (modified) (2 diffs)
• OPA_SRC/TRD/trdmxl.F90 (modified) (2 diffs)
• OPA_SRC/nemogcm.F90 (modified) (3 diffs)
• OPA_SRC/step.F90 (modified) (11 diffs)
• TOP_SRC/TRP/trcadv.F90 (modified) (10 diffs)
• TOP_SRC/TRP/trcldf.F90 (modified) (10 diffs)
• TOP_SRC/TRP/trczdf.F90 (modified) (8 diffs)
• TOP_SRC/trcini.F90 (modified) (6 diffs)

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

@@ -21 +21 @@
 
    IMPLICIT NONE
-   PUBLIC             ! allows the acces to par_oce when dom_oce is used
-   !                  ! exception to coding rules... to be suppressed ???
+   PUBLIC             ! allows the acces to par_oce when dom_oce is used (exception to coding rules)
 
    PUBLIC dom_oce_alloc  ! Called from nemogcm.F90
@@ -107 +106 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdttra  !: vertical profile of tracer time step
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dtra  !: = 2*rdttra except at nit000 (=rdttra) if neuler=0
-
-   !                                         !!* Namelist namcla : cross land advection
-   INTEGER, PUBLIC ::   nn_cla               !: =1 cross land advection for exchanges through some straits (ORCA2)
 
    !!----------------------------------------------------------------------

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

@@ -81 +81 @@
       ENDIF
       !
-                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
+                             CALL dom_nam      ! read namelist ( namrun, namdom )
                              CALL dom_clo      ! Closed seas and lake
                              CALL dom_hgr      ! Horizontal mesh
@@ -131 +131 @@
       !! ** input   : - namrun namelist
       !!              - namdom namelist
-      !!              - namcla namelist
       !!              - namnc4 namelist   ! "key_netcdf4" only
       !!----------------------------------------------------------------------
@@ -146 +145 @@
          &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
          &             ppa2, ppkth2, ppacr2
-      NAMELIST/namcla/ nn_cla
 #if defined key_netcdf4
       NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
@@ -305 +303 @@
       rdth      = rn_rdth
 
-      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
-      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
-905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
-      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
-906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namcla )
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist namcla'
-         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
-      ENDIF
-      IF ( nn_cla .EQ. 1 ) THEN
-         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2 
-            CONTINUE
-         ELSE
-            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
-         ENDIF
-      ENDIF
-
 #if defined key_netcdf4
       !                             ! NetCDF 4 case   ("key_netcdf4" defined)

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

@@ -17 +17 @@
    !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
    !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
+   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask
    !!----------------------------------------------------------------------
 
@@ -129 +130 @@
       !!               tmask_i  : interior ocean mask
       !!----------------------------------------------------------------------
-      !
-      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      INTEGER  ::   ji, jj, jk               ! dummy loop indices
       INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
       INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
@@ -284 +284 @@
       ! 3. Ocean/land mask at wu-, wv- and w points 
       !----------------------------------------------
-      wmask (:,:,1) = tmask(:,:,1) ! ????????
-      wumask(:,:,1) = umask(:,:,1) ! ????????
-      wvmask(:,:,1) = vmask(:,:,1) ! ????????
-      DO jk=2,jpk
-         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
-         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
-         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
+      wmask (:,:,1) = tmask(:,:,1)     ! surface
+      wumask(:,:,1) = umask(:,:,1)
+      wvmask(:,:,1) = vmask(:,:,1)
+      DO jk = 2, jpk                   ! interior values
+         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
+         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
+         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
       END DO
 
@@ -377 +377 @@
       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
          !                                                 ! Increased lateral friction near of some straits
-         IF( nn_cla == 0 ) THEN
-            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
-            ij0 = 101   ;   ij1 = 101
-            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
-            ij0 = 102   ;   ij1 = 102
-            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
-            !
-            !                                ! Bab el Mandeb : partial slip (fmask=1)
-            ij0 =  87   ;   ij1 =  88
-            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
-            ij0 =  88   ;   ij1 =  88
-            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
-            !
-         ENDIF
+         !                                ! Gibraltar strait  : partial slip (fmask=0.5)
+         ij0 = 101   ;   ij1 = 101
+         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
+         ij0 = 102   ;   ij1 = 102
+         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
+         !
+         !                                ! Bab el Mandeb : partial slip (fmask=1)
+         ij0 =  87   ;   ij1 =  88
+         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
+         ij0 =  88   ;   ij1 =  88
+         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
+         !
          !                                ! Danish straits  : strong slip (fmask > 2)
 ! We keep this as an example but it is instable in this case 

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90

@@ -37 +37 @@
       !!                -> not good if located at too high latitude...
       !!----------------------------------------------------------------------
-      !
       REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point
       INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point
@@ -49 +48 @@
       IF( nn_timing == 1 )  CALL timing_start('dom_ngb')
       !
-      CALL wrk_alloc( jpi, jpj, zglam, zgphi, zmask, zdist )
+      CALL wrk_alloc( jpi,jpj,   zglam, zgphi, zmask, zdist )
       !
       zmask(:,:) = 0._wp
@@ -76 +75 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, zglam, zgphi, zmask, zdist )
+      CALL wrk_dealloc( jpi,jpj,   zglam, zgphi, zmask, zdist )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dom_ngb')

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

@@ -36 +36 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-
 
    SUBROUTINE dom_wri_coordinate

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

@@ -501 +501 @@
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
                !                                             ! =====================
-               IF( nn_cla == 0 ) THEN
-                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open 
-                  ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
-                  DO ji = mi0(ii0), mi1(ii1)
-                     DO jj = mj0(ij0), mj1(ij1)
-                        mbathy(ji,jj) = 15
-                     END DO
+               !
+               ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open 
+               ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
+               DO ji = mi0(ii0), mi1(ii1)
+                  DO jj = mj0(ij0), mj1(ij1)
+                     mbathy(ji,jj) = 15
                   END DO
-                  IF(lwp) WRITE(numout,*)
-                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
-                  !
-                  ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
-                  ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
-                  DO ji = mi0(ii0), mi1(ii1)
-                     DO jj = mj0(ij0), mj1(ij1)
-                        mbathy(ji,jj) = 12
-                     END DO
+               END DO
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
+               !
+               ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
+               ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
+               DO ji = mi0(ii0), mi1(ii1)
+                  DO jj = mj0(ij0), mj1(ij1)
+                     mbathy(ji,jj) = 12
                   END DO
-                  IF(lwp) WRITE(numout,*)
-                  IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
-               ENDIF
+               END DO
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
                !
             ENDIF
@@ -545 +544 @@
             !       
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
-               !
-              IF( nn_cla == 0 ) THEN
-                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open 
-                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
-                 DO ji = mi0(ii0), mi1(ii1)
-                    DO jj = mj0(ij0), mj1(ij1)
-                       bathy(ji,jj) = 284._wp
-                    END DO
+            !
+              ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open 
+              ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
+              DO ji = mi0(ii0), mi1(ii1)
+                 DO jj = mj0(ij0), mj1(ij1)
+                    bathy(ji,jj) = 284._wp
                  END DO
-                 IF(lwp) WRITE(numout,*)     
-                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
-                 !
-                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
-                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
-                 DO ji = mi0(ii0), mi1(ii1)
-                    DO jj = mj0(ij0), mj1(ij1)
-                       bathy(ji,jj) = 137._wp
-                    END DO
+               END DO
+              IF(lwp) WRITE(numout,*)     
+              IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
+              !
+              ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
+              ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
+               DO ji = mi0(ii0), mi1(ii1)
+                 DO jj = mj0(ij0), mj1(ij1)
+                    bathy(ji,jj) = 137._wp
                  END DO
-                 IF(lwp) WRITE(numout,*)
-                 IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
-              ENDIF
+              END DO
+              IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
               !
            ENDIF

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

@@ -174 +174 @@
             END DO
          END DO
-         IF( nn_cla == 1 ) THEN                          ! Cross Land advection
-            il0 = 138   ;   il1 = 138                          ! set T & S profile at Gibraltar Strait
-            ij0 = 101   ;   ij1 = 102
-            ii0 = 139   ;   ii1 = 139
-            DO jl = mi0(il0), mi1(il1)
-               DO jj = mj0(ij0), mj1(ij1)
-                  DO ji = mi0(ii0), mi1(ii1)
-                     sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:)
-                     sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:)
-                  END DO
-               END DO
-            END DO
-            il0 = 164   ;   il1 = 164                          ! set T & S profile at Bab el Mandeb Strait
-            ij0 =  87   ;   ij1 =  88
-            ii0 = 161   ;   ii1 = 163
-            DO jl = mi0(il0), mi1(il1)
-               DO jj = mj0(ij0), mj1(ij1)
-                  DO ji = mi0(ii0), mi1(ii1)
-                     sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:)
-                     sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:)
-                  END DO
-               END DO
-            END DO
-         ELSE                                            ! No Cross Land advection
-            ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea
-            ii0 = 148   ;   ii1 = 160
-            sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp
-            sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
-            sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
-         ENDIF
+         ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea
+         ii0 = 148   ;   ii1 = 160
+         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp
+         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
+         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
       ENDIF
       !

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

@@ -34 +34 @@
    USE dtatsd          ! data temperature and salinity   (dta_tsd routine)
    USE dtauvd          ! data: U & V current             (dta_uvd routine)
-   USE zpshde          ! partial step: hor. derivative (zps_hde routine)
-   USE eosbn2          ! equation of state            (eos bn2 routine)
    USE domvvl          ! varying vertical mesh
    USE dynspg_oce      ! pressure gradient schemes
@@ -119 +117 @@
             ENDIF
             IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
-               CALL wrk_alloc( jpi, jpj, jpk, 2, zuvd )
+               CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
                CALL dta_uvd( nit000, zuvd )
                ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
                vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
-               CALL wrk_dealloc( jpi, jpj, jpk, 2, zuvd )
+               CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
             ENDIF
          ENDIF
-         !
+         !   
          ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
          IF( lk_vvl ) THEN
             DO jk = 1, jpk
                fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
-            ENDDO
+            END DO
          ENDIF
          ! 
@@ -456 +454 @@
       !!----------------------------------------------------------------------
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zprn)
+      CALL wrk_alloc( jpi,jpj,jpk,   zprn)
       !
       IF(lwp) WRITE(numout,*) 
@@ -553 +551 @@
       rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zprn)
+      CALL wrk_dealloc( jpi,jpj,jpk,   zprn)
       !
    END SUBROUTINE istate_uvg

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

@@ -18 +18 @@
    !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here
    !!            3.6  ! 2014-11  (P. Mathiot)          isf            added directly here
+   !!            3.7  ! 2015-10  (G. Madec) remove cross-land advection
    !!----------------------------------------------------------------------
 
@@ -29 +30 @@
    USE sbcrnf          ! river runoff 
    USE sbcisf          ! ice shelf 
-   USE cla             ! cross land advection             (cla_div routine)
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -68 +68 @@
       !!         - compute the now divergence given by :
       !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
-      !!      correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf)
-      !!      and cross land flow (div_cla) 
+      !!      correct hdiv with runoff inflow (div_rnf) and ice shelf melting (div_isf)
       !!              II. vorticity :
       !!         - save the curl computed at the previous time-step
@@ -226 +225 @@
       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field)
       IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
-      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
       
       ! 4. Lateral boundary conditions on hdivn and rotn
@@ -256 +254 @@
       !!      - compute the now divergence given by :
       !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
-      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla) 
+      !!      correct hdiv with runoff inflow (div_rnf)
       !!              - Relavtive Vorticity :
       !!      - save the curl computed at the previous time-step (rotb = rotn)
@@ -325 +323 @@
       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field)
       IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
-      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
       !
       CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

@@ -271 +271 @@
       CALL solver_init( nit000 )   ! Elliptic solver initialisation
 #endif
-
-      !                        ! Control of timestep choice
-      IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN
-         IF( nn_cla == 1 )   CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )
-      ENDIF
-
       !               ! Control of hydrostatic pressure choice
-      IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN
-         CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )
-      ENDIF
+      IF( lk_dynspg_ts .AND. ln_dynhpg_imp )   CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init')

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

@@ -14 +14 @@
    !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
    !!            3.7  !  2014-04  (F. Roquet, G. Madec)  add some trends diag
+   !!             -   !  2014-12  (G. Madec) remove cross-land advection (cla)
    !!----------------------------------------------------------------------
 #if defined key_dynspg_flt   ||   defined key_esopa  
@@ -36 +37 @@
    USE bdydyn          ! ocean open boundary condition on dynamics
    USE bdyvol          ! ocean open boundary condition (bdy_vol routine)
-   USE cla             ! cross land advection
    USE trd_oce         ! trends: ocean variables
    USE trddyn          ! trend manager: dynamics
@@ -206 +206 @@
       CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces 
 #endif
-      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
 
       ! compute the next vertically averaged velocity (effect of the additional force not included)

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -80 +80 @@
       IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
       !
-      CALL wrk_alloc( jpi, jpj, zhdiv ) 
+      CALL wrk_alloc( jpi,jpj,   zhdiv ) 
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         !
       ENDIF
       !
@@ -159 +157 @@
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt           ! time step
+      INTEGER, INTENT(in) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   z1_2dt       ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      !
-      INTEGER             ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp)            ::   z1_2dt       ! local scalars
-      !!----------------------------------------------------------------------
-      
+      !!----------------------------------------------------------------------
+      !
       IF( nn_timing == 1 )  CALL timing_start('wzv')
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
@@ -177 +173 @@
          !
          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
       ENDIF
       !                                           !------------------------------!
@@ -216 +211 @@
 
 #if defined key_bdy
-      IF (lk_bdy) THEN
+      IF( lk_bdy ) THEN
          DO jk = 1, jpkm1
             wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
@@ -224 +219 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
-
-
+      !
    END SUBROUTINE wzv
+
 
    SUBROUTINE ssh_swp( kt )
@@ -265 +260 @@
          sshb(:,:) = sshn(:,:)                           ! before <-- now
          sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
+         !
       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

@@ -7 +7 @@
    !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!            3.6  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation
-   !!            4.0  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes 
+   !!            3.7  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes 
    !!             -   !  2014-12  (G. Madec) suppression of cross land advection option
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   tra_adv      : compute ocean tracer advection trend
-   !!   tra_adv_ctl  : control the different options of advection scheme
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers
-   USE dom_oce         ! ocean space and time domain
-   USE domvvl          ! variable vertical scale factors
-   USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine)
-   USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine)
-   USE traadv_muscl    ! MUSCL    scheme           (tra_adv_muscl  routine)
-   USE traadv_muscl2   ! MUSCL2   scheme           (tra_adv_muscl2 routine)
-   USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine)
-   USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine)
-   USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine)
-   USE cla             ! cross land advection      (cla_traadv     routine)
-   USE ldftra          ! lateral diffusion coefficient on tracers
-   USE ldfslp          ! Lateral diffusion: slopes of neutral surfaces
+   !!   tra_adv       : compute ocean tracer advection trend
+   !!   tra_adv_ctl   : control the different options of advection scheme
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE domvvl         ! variable vertical scale factors
+   USE traadv_cen     ! centered scheme           (tra_adv_cen  routine)
+   USE traadv_fct     ! FCT      scheme           (tra_adv_fct  routine)
+   USE traadv_mus     ! MUSCL    scheme           (tra_adv_mus  routine)
+   USE traadv_ubs     ! UBS      scheme           (tra_adv_ubs  routine)
+   USE traadv_qck     ! QUICKEST scheme           (tra_adv_qck  routine)
+   USE traadv_mle     ! ML eddy induced velocity  (tra_adv_mle  routine)
+   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff.
+   USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces
    !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE prtctl          ! Print control
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
-   USE sbc_oce
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O module
+   USE prtctl         ! Print control
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
+
    USE diaptr          ! Poleward heat transport 
-
 
    IMPLICIT NONE
@@ -45 +42 @@
    PUBLIC   tra_adv_init   ! routine called by opa module
 
-   !                              !!* Namelist namtra_adv *
-   LOGICAL ::   ln_traadv_cen2     ! 2nd order centered scheme flag
-   LOGICAL ::   ln_traadv_tvd      ! TVD scheme flag
-   LOGICAL ::   ln_traadv_tvd_zts  ! TVD scheme flag with vertical sub time-stepping
-   LOGICAL ::   ln_traadv_muscl    ! MUSCL scheme flag
-   LOGICAL ::   ln_traadv_muscl2   ! MUSCL2 scheme flag
-   LOGICAL ::   ln_traadv_ubs      ! UBS scheme flag
-   LOGICAL ::   ln_traadv_qck      ! QUICKEST scheme flag
-   LOGICAL ::   ln_traadv_msc_ups  ! use upstream scheme within muscl
-
-
-   INTEGER ::   nadv   ! choice of the type of advection scheme
-
+   !                            !!* Namelist namtra_adv *
+   LOGICAL ::   ln_traadv_cen    ! centered scheme flag
+   INTEGER ::      nn_cen_h, nn_cen_v   ! =2/4 : horizontal and vertical choices of the order of CEN scheme
+   LOGICAL ::   ln_traadv_fct    ! FCT scheme flag
+   INTEGER ::      nn_fct_h, nn_fct_v   ! =2/4 : horizontal and vertical choices of the order of FCT scheme
+   INTEGER ::      nn_fct_zts           ! >=1 : 2nd order FCT with vertical sub-timestepping
+   LOGICAL ::   ln_traadv_mus    ! MUSCL scheme flag
+   LOGICAL ::      ln_mus_ups           ! use upstream scheme in vivcinity of river mouths
+   LOGICAL ::   ln_traadv_ubs    ! UBS scheme flag
+   INTEGER ::      nn_ubs_v             ! =2/4 : vertical choice of the order of UBS scheme
+   LOGICAL ::   ln_traadv_qck    ! QUICKEST scheme flag
+
+   INTEGER ::              nadv             ! choice of the type of advection scheme
+   !
+   !                                        ! associated indices:
+   INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection
+   INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme
+   INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme
+   INTEGER, PARAMETER ::   np_FCT_zts = 3   ! 2nd order FCT scheme with vertical sub-timestepping
+   INTEGER, PARAMETER ::   np_MUS     = 4   ! MUSCL scheme
+   INTEGER, PARAMETER ::   np_UBS     = 5   ! 3rd order Upstream Biased Scheme
+   INTEGER, PARAMETER ::   np_QCK     = 6   ! QUICK scheme
+   
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -84 +91 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_adv')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
+      CALL wrk_alloc( jpi,jpj,jpk,   zun, zvn, zwn )
       !
       !                                          ! set time step
@@ -93 +100 @@
       ENDIF
       !
-      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_traadv( kt )       !==  Cross Land Advection  ==! (hor. advection)
-      !
-      !                                               !==  effective transport  ==!
+      !                                         !==  effective transport  ==!
       DO jk = 1, jpkm1
          zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
@@ -102 +107 @@
       END DO
       !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections
          zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
          zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
       ENDIF
       !
-      zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
-      zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
-      zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
+      zun(:,:,jpk) = 0._wp                                                      ! no transport trough the bottom
+      zvn(:,:,jpk) = 0._wp
+      zwn(:,:,jpk) = 0._wp
       !
       IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   &
@@ -120 +125 @@
       CALL iom_put( "wocetr_eff", zwn )
       !
+!!gm ???
       IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF 
-      !
-   
-      SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
-      CASE ( 1 )   ;    CALL tra_adv_cen2   ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered
-      CASE ( 2 )   ;    CALL tra_adv_tvd    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD 
-      CASE ( 3 )   ;    CALL tra_adv_muscl  ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )   !  MUSCL 
-      CASE ( 4 )   ;    CALL tra_adv_muscl2 ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2 
-      CASE ( 5 )   ;    CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
-      CASE ( 6 )   ;    CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
-      CASE ( 7 )   ;    CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS
-      !
-      CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
-         CALL tra_adv_cen2  ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )          
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv0 - Ta: ', mask1=tmask,               &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_tvd   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv1 - Ta: ', mask1=tmask,               &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )          
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv3 - Ta: ', mask1=tmask,               &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_muscl2( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv4 - Ta: ', mask1=tmask,               &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_ubs   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv5 - Ta: ', mask1=tmask,               &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_qck   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv6 - Ta: ', mask1=tmask,               &
-            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+!!gm ???
+      !
+      SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==!
+      !
+      CASE ( np_CEN )                                    ! Centered scheme : 2nd / 4th order
+         CALL tra_adv_cen    ( kt, nit000, 'TRA',         zun, zvn, zwn     , tsn, tsa, jpts, nn_cen_h, nn_cen_v )
+      CASE ( np_FCT )                                    ! FCT scheme      : 2nd / 4th order
+         CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v )
+      CASE ( np_FCT_zts )                                ! 2nd order FCT with vertical time-splitting
+         CALL tra_adv_fct_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_fct_zts )
+      CASE ( np_MUS )                                    ! MUSCL
+         CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts        , ln_mus_ups ) 
+      CASE ( np_UBS )                                    ! UBS
+         CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts        , nn_ubs_v   )
+      CASE ( np_QCK )                                    ! QUICKEST
+         CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts                     )
+      !
       END SELECT
       !
@@ -159 +152 @@
       IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv' )
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zun, zvn, zwn )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zun, zvn, zwn )
       !                                          
    END SUBROUTINE tra_adv
@@ -171 +164 @@
       !!              tracer advection schemes and set nadv
       !!----------------------------------------------------------------------
-      INTEGER ::   ioptio
-      INTEGER ::   ios                 ! Local integer output status for namelist read
-      !!
-      NAMELIST/namtra_adv/ ln_traadv_cen2 , ln_traadv_tvd,     &
-         &                 ln_traadv_muscl, ln_traadv_muscl2,  &
-         &                 ln_traadv_ubs  , ln_traadv_qck,     &
-         &                 ln_traadv_msc_ups, ln_traadv_tvd_zts
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam_ref )              ! Namelist namtra_adv in reference namelist : Tracer advection scheme
+      INTEGER ::   ioptio, ios   ! Local integers
+      !
+      NAMELIST/namtra_adv/ ln_traadv_cen, nn_cen_h, nn_cen_v,               &   ! CEN
+         &                 ln_traadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts,   &   ! FCT
+         &                 ln_traadv_mus,                     ln_mus_ups,   &   ! MUSCL
+         &                 ln_traadv_ubs,           nn_ubs_v,               &   ! UBS
+         &                 ln_traadv_qck                                        ! QCK
+      !!----------------------------------------------------------------------
+      !
+      !                                !==  Namelist  ==!
+      !
+      REWIND( numnam_ref )                   ! Namelist namtra_adv in reference namelist : Tracer advection scheme
       READ  ( numnam_ref, namtra_adv, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namtra_adv in configuration namelist : Tracer advection scheme
+      !
+      REWIND( numnam_cfg )                   ! Namelist namtra_adv in configuration namelist : Tracer advection scheme
       READ  ( numnam_cfg, namtra_adv, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namtra_adv )
 
-      IF(lwp) THEN                    ! Namelist print
+      IF(lwp) THEN                           ! Namelist print
          WRITE(numout,*)
          WRITE(numout,*) 'tra_adv_init : choice/control of the tracer advection scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_adv : chose a advection scheme for tracers'
-         WRITE(numout,*) '      2nd order advection scheme     ln_traadv_cen2    = ', ln_traadv_cen2
-         WRITE(numout,*) '      TVD advection scheme           ln_traadv_tvd     = ', ln_traadv_tvd
-         WRITE(numout,*) '      MUSCL  advection scheme        ln_traadv_muscl   = ', ln_traadv_muscl
-         WRITE(numout,*) '      MUSCL2 advection scheme        ln_traadv_muscl2  = ', ln_traadv_muscl2
-         WRITE(numout,*) '      UBS    advection scheme        ln_traadv_ubs     = ', ln_traadv_ubs
-         WRITE(numout,*) '      QUICKEST advection scheme      ln_traadv_qck     = ', ln_traadv_qck
-         WRITE(numout,*) '      upstream scheme within muscl   ln_traadv_msc_ups = ', ln_traadv_msc_ups
-         WRITE(numout,*) '      TVD advection scheme with zts  ln_traadv_tvd_zts = ', ln_traadv_tvd_zts
-      ENDIF
-
-      ioptio = 0                      ! Parameter control
-      IF( ln_traadv_cen2   )   ioptio = ioptio + 1
-      IF( ln_traadv_tvd    )   ioptio = ioptio + 1
-      IF( ln_traadv_muscl  )   ioptio = ioptio + 1
-      IF( ln_traadv_muscl2 )   ioptio = ioptio + 1
-      IF( ln_traadv_ubs    )   ioptio = ioptio + 1
-      IF( ln_traadv_qck    )   ioptio = ioptio + 1
-      IF( ln_traadv_tvd_zts)   ioptio = ioptio + 1
-      IF( lk_esopa         )   ioptio =          1
-
-      IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts )   &
-         .AND. ln_isfcav )   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity')
-
-      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' )
-
-      !                              ! Set nadv
-      IF( ln_traadv_cen2   )   nadv =  1
-      IF( ln_traadv_tvd    )   nadv =  2
-      IF( ln_traadv_muscl  )   nadv =  3
-      IF( ln_traadv_muscl2 )   nadv =  4
-      IF( ln_traadv_ubs    )   nadv =  5
-      IF( ln_traadv_qck    )   nadv =  6
-      IF( ln_traadv_tvd_zts)   nadv =  7
-      IF( lk_esopa         )   nadv = -1
-
-      IF(lwp) THEN                   ! Print the choice
+         WRITE(numout,*) '      centered scheme                           ln_traadv_cen = ', ln_traadv_cen
+         WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_fct_h
+         WRITE(numout,*) '            vertical   2nd/4th order               nn_cen_v   = ', nn_fct_v
+         WRITE(numout,*) '      Flux Corrected Transport scheme           ln_traadv_fct = ', ln_traadv_fct
+         WRITE(numout,*) '            horizontal 2nd/4th order               nn_fct_h   = ', nn_fct_h
+         WRITE(numout,*) '            vertical   2nd/4th order               nn_fct_v   = ', nn_fct_v
+         WRITE(numout,*) '            2nd order + vertical sub-timestepping  nn_fct_zts = ', nn_fct_zts
+         WRITE(numout,*) '      MUSCL scheme                              ln_traadv_mus = ', ln_traadv_mus
+         WRITE(numout,*) '            + upstream scheme near river mouths    ln_mus_ups = ', ln_mus_ups
+         WRITE(numout,*) '      UBS scheme                                ln_traadv_ubs = ', ln_traadv_ubs
+         WRITE(numout,*) '            vertical   2nd/4th order               nn_ubs_v   = ', nn_ubs_v
+         WRITE(numout,*) '      QUICKEST scheme                           ln_traadv_qck = ', ln_traadv_qck
+      ENDIF
+
+      ioptio = 0                       !==  Parameter control  ==!
+      IF( ln_traadv_cen )   ioptio = ioptio + 1
+      IF( ln_traadv_fct )   ioptio = ioptio + 1
+      IF( ln_traadv_mus )   ioptio = ioptio + 1
+      IF( ln_traadv_ubs )   ioptio = ioptio + 1
+      IF( ln_traadv_qck )   ioptio = ioptio + 1
+      !
+      IF( ioptio == 0 ) THEN
+         nadv = np_NO_adv
+         CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' )
+      ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' )
+      !
+      IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 )   &          ! Centered
+                        .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 )   ) THEN
+        CALL ctl_stop( 'tra_adv_init: CEN scheme, choose 2nd or 4th order' )
+      ENDIF
+      IF( ln_traadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 )   &          ! FCT
+                        .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 )   ) THEN
+        CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' )
+      ENDIF
+      IF( ln_traadv_fct .AND. nn_fct_zts > 0 ) THEN
+         IF( nn_fct_h == 4 ) THEN
+            nn_fct_h = 2
+            CALL ctl_stop( 'tra_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' )
+         ENDIF
+         IF( lk_vvl ) THEN
+            CALL ctl_stop( 'tra_adv_init: vertical sub-timestepping not allow in non-linear free surface' )
+         ENDIF
+         IF( nn_fct_zts == 1 )   CALL ctl_warn( 'tra_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' )
+      ENDIF
+      IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN     ! UBS
+        CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' )
+      ENDIF
+      IF( ln_traadv_ubs .AND. nn_ubs_v == 4 ) THEN
+         CALL ctl_warn( 'tra_adv_init: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' )
+      ENDIF
+      IF( ln_isfcav ) THEN                                                       ! ice-shelf cavities
+         IF(  ln_traadv_cen .AND. nn_cen_v /= 4    .OR.   &                            ! NO 4th order with ISF
+            & ln_traadv_fct .AND. nn_fct_v /= 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' )
+      ENDIF
+      !
+      !                                !==  used advection scheme  ==!  
+      !                                      ! set nadv
+      IF( ln_traadv_cen                      )   nadv = np_CEN
+      IF( ln_traadv_fct                      )   nadv = np_FCT
+      IF( ln_traadv_fct .AND. nn_fct_zts > 0 )   nadv = np_FCT_zts
+      IF( ln_traadv_mus                      )   nadv = np_MUS
+      IF( ln_traadv_ubs                      )   nadv = np_UBS
+      IF( ln_traadv_qck                      )   nadv = np_QCK
+
+      IF(lwp) THEN                           ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used'
-         IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used'
-         IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used'
-         IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    scheme is used'
-         IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
-         IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
-         IF( nadv ==  7 )   WRITE(numout,*) '         TVD ZTS   scheme is used'
-         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
-      ENDIF
-      !
-      CALL tra_adv_mle_init          ! initialisation of the Mixed Layer Eddy parametrisation (MLE)
+         IF( nadv == np_NO_adv  )   WRITE(numout,*) '         NO T-S advection'
+         IF( nadv == np_CEN     )   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   &
+            &                                                                        ' Vertical   order: ', nn_cen_v
+         IF( nadv == np_FCT     )   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
+            &                                                                        ' Vertical   order: ', nn_fct_v
+         IF( nadv == np_FCT_zts )   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping'
+         IF( nadv == np_MUS     )   WRITE(numout,*) '         MUSCL    scheme is used'
+         IF( nadv == np_UBS     )   WRITE(numout,*) '         UBS      scheme is used'
+         IF( nadv == np_QCK     )   WRITE(numout,*) '         QUICKEST scheme is used'
+      ENDIF
+      !
+      CALL tra_adv_mle_init            !== initialisation of the Mixed Layer Eddy parametrisation (MLE)  ==!
       !
    END SUBROUTINE tra_adv_init

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

@@ -1 @@
-MODULE traadv_tvd
+MODULE traadv_fct
    !!==============================================================================
-   !!                       ***  MODULE  traadv_tvd  ***
-   !! Ocean  tracers:  horizontal & vertical advective trend
+   !!                       ***  MODULE  traadv_fct  ***
+   !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
    !!==============================================================================
-   !! History :  OPA  !  1995-12  (L. Mortier)  Original code
-   !!                 !  2000-01  (H. Loukos)  adapted to ORCA 
-   !!                 !  2000-10  (MA Foujols E.Kestenare)  include file not routine
-   !!                 !  2000-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes
-   !!                 !  2001-07  (E. Durand G. Madec)  adaptation to ORCA config
-   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
-   !!    NEMO    1.0  !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
-   !!            2.0  !  2008-04  (S. Cravatte) add the i-, j- & k- trends computation
-   !!             -   !  2009-11  (V. Garnier) Surface pressure gradient organization
-   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
+   !! History :  3.7  !  2015-09  (L. Debreu, G. Madec)  original code (inspired from traadv_tvd.F90)
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   tra_adv_tvd   : update the tracer trend with the 3D advection trends using a TVD scheme
-   !!   nonosc        : compute monotonic tracer fluxes by a non-oscillatory algorithm 
+   !!  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
@@ -28 +22 @@
    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 in_out_manager ! I/O manager
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tra_adv_tvd        ! routine called by traadv.F90
-   PUBLIC   tra_adv_tvd_zts    ! routine called by traadv.F90
-
-   LOGICAL ::   l_trd   ! flag to compute trends
+   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
@@ -47 +43 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -53 +49 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_tvd ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
-      &                                       ptb, ptn, pta, kjpt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_adv_tvd  ***
+   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  :   TVD scheme, i.e. 2nd order centered scheme with
-      !!       corrected flux (monotonic correction)
-      !!       note: - this advection scheme needs a leap-frog time scheme
+      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
+      !!               (choice through the value of kn_fct)
+      !!               - 4th order compact scheme on the vertical 
+      !!               - corrected flux (monotonic correction) 
       !!
       !! ** Action : - update (pta) with the now advective tracer trends
-      !!             - save the trends 
-      !!----------------------------------------------------------------------
-      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace
-      !
+      !!             - send the trends for further diagnostics
+      !!----------------------------------------------------------------------
       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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
@@ -79 +76 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
       !
-      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
-      INTEGER  ::   ik  
-      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
-      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
-      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz )
+      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices  
+      REAL(wp) ::   z2dtt, 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_tvd : TVD advection scheme on ', cdtype
+         IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-         !
-         l_trd = .FALSE.
-         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
       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.e0   ;    ztrdy(:,:,:) = 0.e0   ;   ztrdz(:,:,:) = 0.e0
+         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
       ENDIF
       !
-      zwi(:,:,:) = 0.e0 ; 
-      !
+      !                                         ! surface & bottom value : flux set to zero one for all
+      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp                ! except at the surface in linear free surface case
+      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp
+      !
+      zwi(:,:,:) = 0._wp        
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         ! 1. Bottom and k=1 value : flux set to zero
-         ! ----------------------------------
-         zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0
-         zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0
-          
-         zwz(:,:,1  ) = 0._wp
-         ! 2. upstream advection with initial mass fluxes & intermediate update
-         ! --------------------------------------------------------------------
-         ! upstream tracer flux in the i and j direction
+         !
+         !        !==  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
@@ -133 +126 @@
             END DO
          END DO
-
-         ! upstream tracer flux in the k direction
-         ! Interior value
-         DO jk = 2, jpkm1
+         !                    !* upstream tracer flux in the k direction *!
+         DO jk = 2, jpkm1         ! Interior value ( multiplied by wmask)
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -145 +136 @@
             END DO
          END DO
-         ! Surface value
-         IF( lk_vvl ) THEN   
-            IF ( ln_isfcav ) THEN
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable
-                  END DO
-               END DO
-            ELSE
-               zwz(:,:,1) = 0.e0          ! volume variable
-            END IF
-         ELSE                
-            IF ( ln_isfcav ) THEN
+         !                    
+         IF(.NOT.lk_vvl ) 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
@@ -163 +144 @@
                   END DO
                END DO   
-            ELSE
-               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
-            END IF
+            ELSE                             ! no cavities: only at the ocean surface
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+            ENDIF
          ENDIF
-
-         ! total advective trend
-         DO jk = 1, jpkm1
+         !               
+         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.
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
                   ! total intermediate advective trends
-                  ztra = - zbtr * (  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) )
+                  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) ) / ( e1e2t(ji,jj) * fse3t_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)
                   zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
@@ -184 +164 @@
             END DO
          END DO
-         !                             ! Lateral boundary conditions on zwi  (unchanged sign)
-         CALL lbc_lnk( zwi, 'T', 1. )  
-
-         !                                 ! trend diagnostics (contribution of upstream fluxes)
-         IF( l_trd )  THEN 
-            ! store intermediate advective trends
+         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)
+         !                             ! "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. antidiffusive flux : high order minus low order
-         ! --------------------------------------------------
-         ! antidiffusive flux on i and j
-         DO jk = 1, jpkm1
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwx(ji,jj,jk) = 0.5 * 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 * 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
-      
-         ! antidiffusive flux on k
-         ! Interior value
-         DO jk = 2, jpkm1                    
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         ! surface value
-         IF ( ln_isfcav ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zwz(ji,jj,mikt(ji,jj)) = 0.e0
-               END DO
-            END DO
-         ELSE
-            zwz(:,:,1) = 0.e0
-         END IF
+         !
+         !
+         !        !==  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                                   ! First 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                                   ! 
+                  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                                 ! gradient
+               DO jj = 1, jpjm1                                   ! First 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
+            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
+         !                                !* vertical anti-diffusive fluxes
+         SELECT CASE( kn_fct_v )                ! Interior values (w-masked)
+         !
+         CASE(  2  )                                  ! 2nd order centered
+            DO jk = 2, jpkm1    
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwz(ji,jj,jk) =  ( 0.5_wp * pwn(ji,jj,jk) * ( 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 )         ! 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
+         !                                      ! top ocean value: high order = upstream  ==>>  zwz=0
+         zwz(:,:, 1 ) = 0._wp                   ! only ocean surface as interior zwz values have been w-masked
+         !
          CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
          CALL lbc_lnk( zwz, 'W',  1. )
 
-         ! 4. monotonicity algorithm
-         ! -------------------------
+         !        !==  monotonicity algorithm  ==!
+         !
          CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
 
 
-         ! 5. final trend with corrected fluxes
-         ! ------------------------------------
+         !        !==  final trend with corrected fluxes  ==!
+         !
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ! total advective trends
-                  ztra = - zbtr * (  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) )
-                  ! add them to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-         !                                 ! trend diagnostics (contribution of upstream fluxes)
-         IF( l_trd )  THEN 
+                  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) ) &
+                     &                                / ( e1e2t(ji,jj) * fse3t(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 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)
@@ -271 +320 @@
       END DO
       !
-                   CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz )
-      IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd')
-      !
-   END SUBROUTINE tra_adv_tvd
-
-
-   SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
-      &                                       ptb, ptn, pta, kjpt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_adv_tvd_zts  ***
+      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 
@@ -296 +344 @@
       !!             - save the trends 
       !!----------------------------------------------------------------------
-      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace
-      !
       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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
@@ -310 +357 @@
       REAL(wp), DIMENSION( jpk )                           ::   zr_p2dt         ! reciprocal of tracer timestep
       INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices  
-      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
       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) ::   z2dtt, zbtr, ztra        ! local scalar
+      REAL(wp) ::   z2dtt, 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, zwz, zhdiv, zwz_sav, zwzts
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd_zts')
-      !
-      CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav )
-      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts )
-      CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs )
+      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_tvd_zts : TVD ZTS advection scheme on ', cdtype
+         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( ( 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 )
+         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( jnzts, wp )
+      z_rzts = 1._wp / REAL( kn_fct_zts, wp )
       zr_p2dt(:) = 1._wp / p2dt(:)
       !
@@ -373 +419 @@
 
          ! upstream tracer flux in the k direction
-         ! Interior value
-         DO jk = 2, jpkm1
+         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) )
-               END DO
-            END DO
-         END DO
-         ! Surface value
-         IF( lk_vvl ) THEN
-            IF ( ln_isfcav ) THEN
+                  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
+         !                       ! top value
+         IF( lk_vvl ) THEN             ! variable volume: only k=1 as zwz is multiplied by wmask
+            zwz(:,:, 1 ) = 0._wp
+         ELSE                          ! linear free surface
+            IF( ln_isfcav ) THEN             ! ice-shelf cavities
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf
-                  END DO
-               END DO
-            ELSE
-               zwz(:,:,1) = 0.e0                              ! volume variable + no isf
-            END IF
-         ELSE
-            IF ( ln_isfcav ) THEN
-               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 +    isf
-                  END DO
-               END DO
-            ELSE
-               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf
-            END IF
+                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 
+                  END DO
+               END DO   
+            ELSE                             ! standard case
+               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+            ENDIF
          ENDIF
-
-         ! total advective trend
-         DO jk = 1, jpkm1
+         !
+         DO jk = 1, jpkm1         ! total advective trend
             z2dtt = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
                   ! total intermediate advective trends
-                  ztra = - zbtr * (  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) )
+                  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) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                   ! update and guess with monotonic sheme
                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra
@@ -422 +457 @@
             END DO
          END DO
-         !                             ! Lateral boundary conditions on zwi  (unchanged sign)
-         CALL lbc_lnk( zwi, 'T', 1. )  
-
-         !                                 ! trend diagnostics (contribution of upstream fluxes)
-         IF( l_trd )  THEN 
-            ! store intermediate advective trends
+         !                           
+         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)
+         !                                ! "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(:,:,:) )
@@ -436 +469 @@
          ENDIF
 
-         ! 3. antidiffusive flux : high order minus low order
-         ! --------------------------------------------------
-         ! antidiffusive flux on i and j
-
-
-         DO jk = 1, jpkm1
-
+         ! 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 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)   &
@@ -459 +490 @@
                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)
+                  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
       
-         ! antidiffusive flux on k
-         zwz(:,:,1) = 0._wp        ! Surface value
-         zwz_sav(:,:,:) = zwz(:,:,:)
-         !
-         ztrs(:,:,:,1) = ptb(:,:,:,jn)
-         zwzts(:,:,:) = 0._wp
-
-         DO jl = 1, jnzts                   ! 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( jnzts + 1 , 2)  ! Toggle to collect every second flux
-                                            ! starting at jl =1 if jnzts 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
+         !                                !* vertical anti-diffusive flux
+         zwz_sav(:,:,:)   = zwz(:,:,:)
+         ztrs   (:,:,:,1) = ptb(:,:,:,jn)
+         zwzts  (:,:,:)   = 0._wp
+         IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp    ! surface value set to zero in vvl case
+         !
+         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 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) )
-                     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
-
+                     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(.NOT.lk_vvl ) 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                                      ! standard case
+                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+               ENDIF
+            ENDIF
+            !
             jtaken = MOD( jtaken + 1 , 2 )
-
-            DO jk = 2, jpkm1          ! Interior value
+            !
+            DO jk = 2, jpkm1                             ! total advective trends
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
-                     ! total advective trends
-                     ztra = - zbtr * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
-                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra
-                  END DO
-               END DO
-            END DO
-
+                     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) )   &
+                        &                         / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
+                  END DO
+               END DO
+            END DO
+            !
          END DO
 
@@ -518 +558 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk)
+                  zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)
                END DO
             END DO
@@ -535 +575 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ! total advective trends
-                  ztra = - zbtr * (  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) )
-                  ! add them to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  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)   )   &
+                     &                                / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                END DO
             END DO
@@ -551 +587 @@
             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)
@@ -564 +602 @@
       END DO
       !
-                   CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts )
-                   CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs )
-                   CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav )
-      IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd_zts')
-      !
-   END SUBROUTINE tra_adv_tvd_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( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts')
+      !
+   END SUBROUTINE tra_adv_fct_zts
 
 
@@ -683 +720 @@
    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_tvd
+END MODULE traadv_fct

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90

@@ -1 @@
-MODULE traadv_muscl
+MODULE traadv_mus
    !!======================================================================
-   !!                       ***  MODULE  traadv_muscl  ***
+   !!                       ***  MODULE  traadv_mus  ***
    !! Ocean  tracers:  horizontal & vertical advective trend
    !!======================================================================
@@ -9 +9 @@
    !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
    !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   tra_adv_muscl : update the tracer trend with the horizontal
+   !!            3.7  !  2015-09  (G. Madec) add the ice-shelf cavities boundary condition
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_adv_mus   : update the tracer trend with the horizontal
    !!                   and vertical advection trends using MUSCL scheme
    !!----------------------------------------------------------------------
@@ -34 +35 @@
    PRIVATE
 
-   PUBLIC   tra_adv_muscl   ! routine called by traadv.F90
+   PUBLIC   tra_adv_mus   ! routine called by traadv.F90
    
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
@@ -50 +51 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn,   &
-      &                                                ptb, pta, kjpt, ld_msc_ups )
+   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn,             &
+      &                                              ptb, pta, kjpt, ld_msc_ups )
       !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE tra_adv_muscl  ***
-      !!
-      !! ** Purpose :   Compute the now trend due to total advection of T and 
-      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
-      !!      Conservation Laws) and add it to the general tracer trend.
+      !!                    ***  ROUTINE tra_adv_mus  ***
+      !!
+      !! ** Purpose :   Compute the now trend due to total advection of tracers
+      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for
+      !!              Conservation Laws) and add it to the general tracer trend.
       !!
       !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
+      !!              ld_msc_ups=T : 
       !!
       !! ** Action  : - update (ta,sa) with the now advective tracer trends
-      !!              - save trends 
+      !!              - send trends to trdtra module for further diagnostcs
       !!
       !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
@@ -77 +79 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend 
       !
-      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
-      INTEGER  ::   ierr                      ! local integer
-      REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars
-      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      -
-      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
+      INTEGER  ::   ji, jj, jk, jn       ! dummy loop indices
+      INTEGER  ::   ierr                 ! local integer
+      REAL(wp) ::   zu, z0u, zzwx, zw    ! local scalars
+      REAL(wp) ::   zv, z0v, zzwy, z0w   !   -      -
+      REAL(wp) ::   zdt, zalpha          !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      - 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_adv_muscl')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy )
+      IF( nn_timing == 1 )  CALL timing_start('tra_adv_mus')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy )
       !
       IF( kt == kit000 )  THEN
@@ -97 +99 @@
          IF(lwp) WRITE(numout,*)
          !
-         !
-         IF( ld_msc_ups ) THEN
-            IF( .NOT. ALLOCATED( upsmsk ) )  THEN
-                ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
-                IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate upsmsk array')
-            ENDIF
+         ! Upstream / MUSCL scheme indicator
+         !
+         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
+         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
+         !
+         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
+            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
             upsmsk(:,:) = 0._wp                             ! not upstream by default
-         ENDIF
-
-         IF( .NOT. ALLOCATED( xind ) ) THEN
-             ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
-             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array')
-         ENDIF
-         !
-
-         !
-         ! Upstream / MUSCL scheme indicator
-         ! ------------------------------------
-!!gm  useless
-         xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed
-!!gm
-         !
-         IF( ld_msc_ups )  THEN
+            !
             DO jk = 1, jpkm1
                xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
                   &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
-                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 near some straits
+                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
             END DO
          ENDIF 
@@ -172 +160 @@
                END DO
            END DO
-         END DO             ! interior values
-
+         END DO
+         !
          !                                             !-- MUSCL horizontal advective fluxes
          DO jk = 1, jpkm1                                     ! interior values
@@ -196 +184 @@
             END DO
          END DO
-         !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign)
-         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
-         !
-         ! Tracer flux divergence at t-point added to the general trend
-         DO jk = 1, jpkm1
+         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
+         !
+         DO jk = 1, jpkm1         ! Tracer flux divergence at t-point added to the general trend
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ! horizontal advective trends
-                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                  &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
-                  ! add it to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  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  ) )     &
+                  &                                   / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
                END DO
            END DO
@@ -226 +209 @@
          ! II. Vertical advective fluxes
          ! -----------------------------
-         !                                             !-- first guess of the slopes
-         zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions
+         !                                !-- first guess of the slopes
+         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
+         zwx(:,:,jpk) = 0._wp                   ! surface & bottom boundary conditions
          DO jk = 2, jpkm1                                     ! interior values
             zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
          END DO
 
-         !                                             !-- Slopes of tracer
-         zslpx(:,:,1) = 0.e0                                  ! surface values
-         DO jk = 2, jpkm1                                     ! interior value
+         !                                !-- Slopes of tracer
+         zslpx(:,:,1) = 0._wp                   ! surface values
+         DO jk = 2, jpkm1                       ! interior value
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -242 +226 @@
             END DO
          END DO
-         !                                             !-- Slopes limitation
-         DO jk = 2, jpkm1                                     ! interior values
+         !                                !-- Slopes limitation
+         DO jk = 2, jpkm1                       ! interior values
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -252 +236 @@
             END DO
          END DO
-         !                                             !-- vertical advective flux
-         !                                                    ! surface values  (bottom already set to zero)
-         IF( lk_vvl ) THEN    ;   zwx(:,:, 1 ) = 0.e0                      !  variable volume
-         ELSE                 ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
-         ENDIF 
-         !
-         DO jk = 1, jpkm1                                     ! interior values
+         !                                !-- vertical advective flux
+         DO jk = 1, jpkm1                       ! interior values
             zdt  = p2dt(jk)
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) )
                   z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
                   zalpha = 0.5 + z0w
-                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 
+                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt  / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) )
                   zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
                   zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
-                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
                END DO 
             END DO
          END DO
-
+         !                                      ! top values  (bottom already set to zero)
+         IF( lk_vvl ) THEN                            !  variable volume
+            zwx(:,:, 1 ) = 0._wp                            ! k=1 only as zwx has been multiplied by wmask
+         ELSE                                         ! linear free surface 
+            IF( ln_isfcav ) THEN                            ! ice-shelf cavities (top of the ocean)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwx(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: only at the ocean surface
+               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
+            ENDIF
+         ENDIF
+         !
          DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ! vertical advective trends 
-                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
-                  ! add it to the general tracer trends
-                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) + ztra
+                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
                END DO
             END DO
@@ -291 +279 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl')
-      !
-   END SUBROUTINE tra_adv_muscl
+      CALL wrk_dealloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mus')
+      !
+   END SUBROUTINE tra_adv_mus
 
    !!======================================================================
-END MODULE traadv_muscl
+END MODULE traadv_mus

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

@@ -102 +102 @@
          IF(lwp) WRITE(numout,*)
       ENDIF
+      !
       l_trd = .FALSE.
       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
@@ -130 +131 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
       !!
-      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zfu, zfc, zfd
       !----------------------------------------------------------------------
       !
@@ -139 +140 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0  
-         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
-         !                                                  
-         DO jk = 1, jpkm1                                
-            !                                             
-            !--- Computation of the ustream and downstream value of the tracer and the mask
+         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
+         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
+         !
+!!gm why not using a SHIFT instruction...
+         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! Upstream in the x-direction for the tracer
-                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
-                  ! Downstream in the x-direction for the tracer
-                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
+                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
+                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
                END DO
             END DO
@@ -159 +157 @@
          ! Horizontal advective fluxes
          ! ---------------------------
-         !
          DO jk = 1, jpkm1                             
             DO jj = 2, jpjm1
@@ -380 +377 @@
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   zbtr , ztra      ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
       !!----------------------------------------------------------------------
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zwz )
+      CALL wrk_alloc( jpi,jpj,jpk,   zwz )
+      !
+      !                          ! surface & bottom values 
+      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all
+                     zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface
+      !
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-         ! 1. Bottom value : flux set to zero
-         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
-         !
-         !                                 ! Surface value
-         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
-         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
+         !
+         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         IF(.NOT.lk_vvl ) THEN               !* top value   (only in linear free surf. as zwz is multiplied by wmask)
+            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface 
+                  END DO
+               END DO   
+            ELSE                                   ! no ice-shelf cavities (only ocean surface)
+               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
+            ENDIF
          ENDIF
          !
-         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
+         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
-               END DO
-            END DO
-         END DO
-         !
-         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  ! k- vertical advective trends 
-                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
-                  ! added to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
+                     &                                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
                END DO
             END DO

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

@@ -16 +16 @@
    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
@@ -38 +39 @@
 #  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 +45 @@
 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 +53 @@
       !!      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 +62 @@
       !!      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 +85 @@
       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 +100 @@
       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 +111 @@
       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 = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
@@ -124 +128 @@
                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 +136 @@
             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 +167 @@
             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 +181 @@
             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._wp / ( e1e2t(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 +274 @@
       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 +307 @@
       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 +329 @@
          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 +343 @@
          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 +357 @@
                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 +363 @@
          END DO
       END DO
+      !
       ! monotonic flux in the k direction, i.e. pcc
       ! -------------------------------------------
@@ -366 +377 @@
       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')

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

@@ -73 +73 @@
       SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend
       !
-      CASE ( n_lap   )                                   ! laplacian: iso-level operator
+      CASE ( np_lap   )                                  ! laplacian: iso-level operator
          CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   )
          !
-      CASE ( n_lap_i )                                   ! laplacian: standard iso-neutral operator (Madec)
+      CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec)
          CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
          !
-      CASE ( n_lap_it )                                  ! laplacian: triad iso-neutral operator (griffies)
+      CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies)
          CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
          !
-      CASE ( n_blp , n_blp_i , n_blp_it )                ! bilaplacian: iso-level & iso-neutral operators
+      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
          CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf )
       END SELECT
@@ -126 +126 @@
       IF( ln_traldf_blp )   ioptio = ioptio + 1
       IF( ioptio >  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on tracer' )
-      IF( ioptio == 0   )   nldf = n_no_ldf     ! No lateral diffusion
+      IF( ioptio == 0   )   nldf = np_no_ldf     ! No lateral diffusion
       ioptio = 0
       IF( ln_traldf_lev )   ioptio = ioptio + 1
@@ -137 +137 @@
       IF( ln_traldf_lap ) THEN         ! laplacian operator
          IF ( ln_zco ) THEN               ! z-coordinate
-            IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level = horizontal (no rotation)
-            IF ( ln_traldf_hor   )   nldf = n_lap     ! iso-level = horizontal (no rotation)
-            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard  (   rotation)
-            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad     (   rotation)
+            IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_hor   )   nldf = np_lap     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation)
+            IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation)
          ENDIF
          IF ( ln_zps ) THEN               ! z-coordinate with partial step
-            IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed 
-            IF ( ln_traldf_hor   )   nldf = n_lap     ! horizontal (no rotation)
-            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation)
-            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_traldf_lev   )   ierr = 1          ! iso-level not allowed 
+            IF ( ln_traldf_hor   )   nldf = np_lap     ! horizontal             (no rotation)
+            IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard     (rotation)
+            IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad        (rotation)
          ENDIF
          IF ( ln_sco ) THEN               ! s-coordinate
-            IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level  (no rotation)
-            IF ( ln_traldf_hor   )   nldf = n_lap_i   ! horizontal (   rotation)   
-            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation)
-            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_traldf_lev   )   nldf = np_lap     ! iso-level              (no rotation)
+            IF ( ln_traldf_hor   )   nldf = np_lap_i   ! horizontal             (   rotation)
+            IF ( ln_traldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation)
+            IF ( ln_traldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation)
          ENDIF
       ENDIF
@@ -158 +158 @@
       IF( ln_traldf_blp ) THEN         ! bilaplacian operator
          IF ( ln_zco ) THEN               ! z-coordinate
-            IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level = horizontal (no rotation)
-            IF ( ln_traldf_hor   )   nldf = n_blp     ! iso-level = horizontal (no rotation)
-            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
-            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_hor   )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation)
+            IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation)
          ENDIF
          IF ( ln_zps ) THEN               ! z-coordinate with partial step
-            IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed 
-            IF ( ln_traldf_hor   )   nldf = n_blp     ! horizontal (no rotation)
-            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
-            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_traldf_lev   )   ierr = 1          ! iso-level not allowed 
+            IF ( ln_traldf_hor   )   nldf = np_blp     ! horizontal             (no rotation)
+            IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation)
+            IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation)
          ENDIF
          IF ( ln_sco ) THEN               ! s-coordinate
-            IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level  (no rotation)
-            IF ( ln_traldf_hor   )   nldf = n_blp_it  ! horizontal (   rotation)       !!gm   a checker....
-            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
-            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_traldf_lev   )   nldf = np_blp     ! iso-level              (no rotation)
+            IF ( ln_traldf_hor   )   nldf = np_blp_it  ! horizontal             (   rotation)
+            IF ( ln_traldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard  (   rotation)
+            IF ( ln_traldf_triad )   nldf = np_blp_it  ! iso-neutral: triad     (   rotation)
          ENDIF
       ENDIF
@@ -181 +181 @@
            &            CALL ctl_stop( '          eddy induced velocity on tracers requires isopycnal',    &
            &                                                                    ' laplacian diffusion' )
-      IF(  nldf == n_lap_i .OR. nldf == n_lap_it .OR. &
-         & nldf == n_blp_i .OR. nldf == n_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required 
+      IF(  nldf == np_lap_i .OR. nldf == np_lap_it .OR. &
+         & nldf == np_blp_i .OR. nldf == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required 
       !
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nldf == n_no_ldf )   WRITE(numout,*) '          NO lateral diffusion'
-         IF( nldf == n_lap    )   WRITE(numout,*) '          laplacian iso-level operator'
-         IF( nldf == n_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)'
-         IF( nldf == n_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)'
-         IF( nldf == n_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator'
-         IF( nldf == n_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)'
-         IF( nldf == n_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)'
+         IF( nldf == np_no_ldf )   WRITE(numout,*) '          NO lateral diffusion'
+         IF( nldf == np_lap    )   WRITE(numout,*) '          laplacian iso-level operator'
+         IF( nldf == np_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)'
+         IF( nldf == np_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)'
+         IF( nldf == np_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator'
+         IF( nldf == np_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)'
+         IF( nldf == np_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)'
       ENDIF
       !

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_blp.F90

@@ -8 +8 @@
 
    !!----------------------------------------------------------------------
-   !!   traldf_blp   : update the tracer trend with the bilaplacian lateral mixing trend
+   !!   traldf_blp    : update the tracer trend with the bilaplacian lateral mixing trend
    !!----------------------------------------------------------------------
-   USE oce           ! ocean dynamics and active tracers
-   USE dom_oce       ! ocean space and time domain
-   USE phycst        ! physical constants
-   USE trc_oce       ! share passive tracers/Ocean variables
-   USE zdf_oce       ! ocean vertical physics
-   USE ldftra        ! lateral physics: eddy diffusivity
-   USE ldfslp        ! lateral physics: iso-neutral slopes
-   USE traldf_lap    ! lateral diffusion (Standard operator)         (tra_ldf_lap   routine)
-   USE traldf_iso    ! lateral diffusion (Standard operator)         (tra_ldf_iso   routine)
-   USE traldf_triad  ! lateral diffusion (Standard operator)         (tra_ldf_triad routine)
-   USE diaptr        ! poleward transport diagnostics
-   USE zpshde        ! partial step: hor. derivative     (zps_hde routine)
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE trc_oce        ! share passive tracers/Ocean variables
+   USE zdf_oce        ! ocean vertical physics
+   USE ldftra         ! lateral physics: eddy diffusivity
+   USE ldfslp         ! lateral physics: iso-neutral slopes
+   USE traldf_lap     ! lateral diffusion (Standard operator)         (tra_ldf_lap   routine)
+   USE traldf_iso     ! lateral diffusion (Standard operator)         (tra_ldf_iso   routine)
+   USE traldf_triad   ! lateral diffusion (Standard operator)         (tra_ldf_triad routine)
+   USE diaptr         ! poleward transport diagnostics
+   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
    !
    USE in_out_manager ! I/O manager
@@ -33 +33 @@
    PRIVATE
 
-   PUBLIC   tra_ldf_blp         ! routine called by traldf.F90
+   PUBLIC   tra_ldf_blp   ! routine called by traldf.F90
 
-   !                                    ! Flag to control the type of lateral diffusive operator
-   INTEGER, PARAMETER, PUBLIC ::   n_no_ldf = 00                      ! without operator (i.e. no lateral diffusive trend)
-   !               !!      laplacian    !    bilaplacian   !
-   INTEGER, PARAMETER, PUBLIC ::   n_lap    = 10   ,   n_blp    = 20  ! iso-level operator
-   INTEGER, PARAMETER, PUBLIC ::   n_lap_i  = 11   ,   n_blp_i  = 21  ! standard iso-neutral or geopotential operator 
-   INTEGER, PARAMETER, PUBLIC ::   n_lap_it = 12   ,   n_blp_it = 22  ! triad    iso-neutral or geopotential operator
+   !                      ! Flag to control the type of lateral diffusive operator
+   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral diffusive trend)
+   !                          !!      laplacian     !    bilaplacian    !
+   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator
+   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator
+   INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator
 
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels
@@ -94 +94 @@
          WRITE(numout,*)
          SELECT CASE ( kldf )
-         CASE ( n_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype
-         CASE ( n_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
-         CASE ( n_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
+         CASE ( np_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype
+         CASE ( np_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
+         CASE ( np_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
          END SELECT
          WRITE(numout,*) '~~~~~~~~~~~'
@@ -105 +105 @@
       SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==!
       !
-      CASE ( n_blp    )                ! iso-level bilaplacian
+      CASE ( np_blp    )               ! iso-level bilaplacian
          CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 )
          !
-      CASE ( n_blp_i  )                ! rotated   bilaplacian : standard operator (Madec)
+      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
          !
-      CASE ( n_blp_it )                ! rotated  bilaplacian : triad operator (griffies)
+      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
          !
@@ -126 +126 @@
       SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==!
       !
-      CASE ( n_blp    )                ! iso-level bilaplacian
+      CASE ( np_blp    )               ! iso-level bilaplacian
          CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 )
          !
-      CASE ( n_blp_i  )                ! rotated   bilaplacian : standard operator (Madec)
+      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
          !
-      CASE ( n_blp_it )                ! rotated  bilaplacian : triad operator (griffies)
+      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
          !

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

@@ -83 +83 @@
       CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
       CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            tsb, tsa, jpts )  !   implicit scheme 
-      CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf0 - Ta: ', mask1=tmask,               &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            tsb, tsa, jpts ) 
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf1 - Ta: ', mask1=tmask,               &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       END SELECT
 !!gm WHY here !   and I don't like that !
@@ -149 +142 @@
             &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' )
 
-      ! Test: esopa
-      IF( lk_esopa )    nzdf = -1                      ! All schemes used
-
       IF(lwp) THEN
          WRITE(numout,*)
          WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
          IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
          IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90

@@ -800 +800 @@
       END IF
 
-      IF( nn_cla == 1 )   CALL ctl_warn( '      You set n_cla = 1. Note that the Mixed-Layer diagnostics  ',   &
-         &                               '      are not exact along the corresponding straits.            ')
-
       !                                   ! allocate trdmxl arrays
       IF( trd_mxl_alloc()    /= 0 )   CALL ctl_stop( 'STOP', 'trd_mxl_init : unable to allocate trdmxl     arrays' )
@@ -809 +806 @@
 
 
-
-       nkstp     = nit000 - 1              ! current time step indicator initialization
+      nkstp     = nit000 - 1              ! current time step indicator initialization
 
 

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -33 +33 @@
    !!             -   ! 2013-06  (I. Epicoco, S. Mocavero, CMCC) nemo_northcomms: setup avoiding MPI communication 
    !!            3.7  ! 2014-12  (G. Madec) suppression of cross land advection option
-   !!             -   ! 2014-12  (G. Madec) remove KPP scheme
+   !!             -   ! 2014-12  (G. Madec) remove KPP scheme and cross-land advection (cla)
    !!----------------------------------------------------------------------
 
@@ -46 +46 @@
    !!----------------------------------------------------------------------
    USE step_oce        ! module used in the ocean time stepping module
-   USE cla             ! cross land advection               (tra_cla routine)
    USE domcfg          ! domain configuration               (dom_cfg routine)
    USE mppini          ! shared/distributed memory setting (mpp_init routine)
@@ -458 +457 @@
 
       !                                      ! Misc. options
-      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_init       ! Cross Land Advection
                             CALL sto_par_init   ! Stochastic parametrization
       IF( ln_sto_eos     )  CALL sto_pts_init   ! RRandom T/S fluctuations

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -109 +109 @@
       IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
                          CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
-CALL FLUSH    ( numout )
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
                          CALL sto_par( kstp )          ! Stochastic parameters
-CALL FLUSH    ( numout )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -124 +123 @@
                          CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
                          CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
-CALL FLUSH    ( numout )
+
       !
       !  VERTICAL PHYSICS
@@ -137 +136 @@
          avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
       ENDIF
-CALL FLUSH    ( numout )
+
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
          DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
@@ -152 +151 @@
       IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
       IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
-CALL FLUSH    ( numout )
       !
       !  LATERAL  PHYSICS
@@ -179 +177 @@
       !                                               ! eddy diffusivity coeff. and/or eiv coeff.
       IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp ) 
-write(*,*) 'after ldf_slp'
-CALL FLUSH    ( numout )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -188 +184 @@
       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
                          CALL wzv           ( kstp )  ! now cross-level velocity 
-write(*,*) 'after wzv'
-CALL FLUSH    ( numout )
 
       IF( lk_dynspg_ts ) THEN 
@@ -231 +225 @@
                                   CALL wzv           ( kstp )  ! now cross-level velocity 
       ENDIF
-write(*,*) 'after wzv 2'
-CALL FLUSH    ( numout )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -253 +245 @@
                          CALL trc_stp( kstp )         ! time-stepping
 #endif
-write(*,*) 'end dyn '
-CALL FLUSH    ( numout )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -270 +260 @@
       IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
-write(*,*) 'before tra_ldf'
-CALL FLUSH    ( numout )
                              CALL tra_ldf    ( kstp )       ! lateral mixing
-write(*,*) 'after tra_ldf'
-CALL FLUSH    ( numout )
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
@@ -316 +302 @@
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
       ENDIF
-write(*,*) 'after tra_nxt'
-CALL FLUSH    ( numout )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -360 +344 @@
       IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
       !
-write(*,*) 'after dom_vvl'
-CALL FLUSH    ( numout )
-
 
 !!gm : This does not only concern the dynamics ==>>> add a new title

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90

@@ -4 +4 @@
    !! Ocean passive tracers:  advection trend 
    !!==============================================================================
-   !! History :  2.0  !  05-11  (G. Madec)  Original code
-   !!            3.0  !  10-06  (C. Ethe)   Adapted to passive tracers
+   !! History :  2.0  !  2005-11  (G. Madec)  Original code
+   !!            3.0  !  2010-06  (C. Ethe)   Adapted to passive tracers
+   !!            3.7  !  2014-05  (G. Madec, C. Ethe)  Add 2nd/4th order cases for CEN and FCT schemes 
    !!----------------------------------------------------------------------
 #if defined key_top
@@ -11 +12 @@
    !!   'key_top'                                                TOP models
    !!----------------------------------------------------------------------
-   !!   trc_adv      : compute ocean tracer advection trend
-   !!   trc_adv_ini  : control the different options of advection scheme
-   !!----------------------------------------------------------------------
-   USE oce_trc         ! ocean dynamics and active tracers
-   USE trc             ! ocean passive tracers variables
-   USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine)
-   USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine)
-   USE traadv_muscl    ! MUSCL    scheme           (tra_adv_muscl  routine)
-   USE traadv_muscl2   ! MUSCL2   scheme           (tra_adv_muscl2 routine)
-   USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine)
-   USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine)
-   USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine)
-   USE ldftra          ! lateral diffusion coefficient on tracers
-   USE prtctl_trc      ! Print control
+   !!   trc_adv       : compute ocean tracer advection trend
+   !!   trc_adv_ini   : control the different options of advection scheme
+   !!----------------------------------------------------------------------
+   USE oce_trc        ! ocean dynamics and active tracers
+   USE trc            ! ocean passive tracers variables
+   USE traadv_cen     ! centered scheme           (tra_adv_cen  routine)
+   USE traadv_fct     ! FCT      scheme           (tra_adv_fct  routine)
+   USE traadv_mus     ! MUSCL    scheme           (tra_adv_mus  routine)
+   USE traadv_ubs     ! UBS      scheme           (tra_adv_ubs  routine)
+   USE traadv_qck     ! QUICKEST scheme           (tra_adv_qck  routine)
+   USE traadv_mle     ! ML eddy induced velocity  (tra_adv_mle  routine)
+   USE ldftra         ! lateral diffusion coefficient on tracers
+   USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces
+   !
+   USE prtctl_trc     ! Print control
 
    IMPLICIT NONE
@@ -33 +35 @@
    PUBLIC   trc_adv_ini  
 
-   !                                        !!: ** Advection (namtrc_adv) **
-   LOGICAL , PUBLIC ::   ln_trcadv_cen2      ! 2nd order centered scheme flag
-   LOGICAL , PUBLIC ::   ln_trcadv_tvd       ! TVD scheme flag
-   LOGICAL , PUBLIC ::   ln_trcadv_muscl     ! MUSCL scheme flag
-   LOGICAL , PUBLIC ::   ln_trcadv_muscl2    ! MUSCL2 scheme flag
-   LOGICAL , PUBLIC ::   ln_trcadv_ubs       ! UBS scheme flag
-   LOGICAL , PUBLIC ::   ln_trcadv_qck       ! QUICKEST scheme flag
-   LOGICAL , PUBLIC ::   ln_trcadv_msc_ups   ! use upstream scheme within muscl
-
-   INTEGER ::   nadv   ! choice of the type of advection scheme
+   !                            !!* Namelist namtrc_adv *
+   LOGICAL ::   ln_trcadv_cen    ! centered scheme flag
+   INTEGER ::      nn_cen_h, nn_cen_v   ! =2/4 : horizontal and vertical choices of the order of CEN scheme
+   LOGICAL ::   ln_trcadv_fct    ! FCT scheme flag
+   INTEGER ::      nn_fct_h, nn_fct_v   ! =2/4 : horizontal and vertical choices of the order of FCT scheme
+   INTEGER ::      nn_fct_zts           ! >=1 : 2nd order FCT with vertical sub-timestepping
+   LOGICAL ::   ln_trcadv_mus    ! MUSCL scheme flag
+   LOGICAL ::      ln_mus_ups           ! use upstream scheme in vivcinity of river mouths
+   LOGICAL ::   ln_trcadv_ubs    ! UBS scheme flag
+   INTEGER ::      nn_ubs_v             ! =2/4 : vertical choice of the order of UBS scheme
+   LOGICAL ::   ln_trcadv_qck    ! QUICKEST scheme flag
+
+   !                                        ! choices of advection scheme:
+   INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection
+   INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme
+   INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme
+   INTEGER, PARAMETER ::   np_FCT_zts = 3   ! 2nd order FCT scheme with vertical sub-timestepping
+   INTEGER, PARAMETER ::   np_MUS     = 4   ! MUSCL scheme
+   INTEGER, PARAMETER ::   np_UBS     = 5   ! 3rd order Upstream Biased Scheme
+   INTEGER, PARAMETER ::   np_QCK     = 6   ! QUICK scheme
+
+   INTEGER ::              nadv             ! chosen advection scheme
+   !
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dt  ! vertical profile time-step, = 2 rdttra
    !                                                    ! except at nitrrc000 (=rdttra) if neuler=0
@@ -50 +65 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! NEMO/TOP 3.7 , NEMO Consortium (2015)
    !! $Id$ 
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -60 +75 @@
       !!                  ***  ROUTINE trc_adv_alloc  ***
       !!----------------------------------------------------------------------
-
+      !
       ALLOCATE( r2dt(jpk), STAT=trc_adv_alloc )
-
+      !
       IF( trc_adv_alloc /= 0 ) CALL ctl_warn('trc_adv_alloc : failed to allocate array.')
-
+      !
    END FUNCTION trc_adv_alloc
 
@@ -78 +93 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER ::   jk 
+      INTEGER ::   jk   ! dummy loop index
       CHARACTER (len=22) ::   charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn  ! effective velocity
@@ -93 +108 @@
       ENDIF
       !                                               !==  effective transport  ==!
-      zun(:,:,jpk) = 0._wp                                                       ! no transport trough the bottom
-      zvn(:,:,jpk) = 0._wp
-      zwn(:,:,jpk) = 0._wp
       DO jk = 1, jpkm1
          zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                   ! eulerian transport
@@ -112 +124 @@
       IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )  ! add the mle transport
       !
-      !
-      SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
-      !
-      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
-      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD 
-      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra, ln_trcadv_msc_ups )   !  MUSCL 
-      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2 
-      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS 
-      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST 
+      zun(:,:,jpk) = 0._wp                                                       ! no transport trough the bottom
+      zvn(:,:,jpk) = 0._wp
+      zwn(:,:,jpk) = 0._wp
+      !
+      !
+      SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==!
+      !
+      CASE ( np_CEN )                                    ! Centered : 2nd / 4th order
+         CALL tra_adv_cen    ( kt, nittrc000,'TRC',       zun, zvn, zwn     , trn, tra, jptra, nn_cen_h, nn_cen_v )
+      CASE ( np_FCT )                                    ! FCT      : 2nd / 4th order
+         CALL tra_adv_fct    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra, nn_fct_h, nn_fct_v )
+      CASE ( np_FCT_zts )                                ! 2nd order FCT with vertical time-splitting
+         CALL tra_adv_fct_zts( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra        , nn_fct_zts )
+      CASE ( np_MUS )                                    ! MUSCL
+         CALL tra_adv_mus    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra        , ln_mus_ups ) 
+      CASE ( np_UBS )                                    ! UBS
+         CALL tra_adv_ubs    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra        , nn_ubs_v   )
+      CASE ( np_QCK )                                    ! QUICKEST
+         CALL tra_adv_qck    ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra                     )
       !
       END SELECT
@@ -146 +168 @@
       INTEGER ::  ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namtrc_adv/ ln_trcadv_cen2 , ln_trcadv_tvd   ,    &
-         &                 ln_trcadv_muscl, ln_trcadv_muscl2,    &
-         &                 ln_trcadv_ubs  , ln_trcadv_qck, ln_trcadv_msc_ups
+      NAMELIST/namtrc_adv/ ln_trcadv_cen, nn_cen_h, nn_cen_v,               &   ! CEN
+         &                 ln_trcadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts,   &   ! FCT
+         &                 ln_trcadv_mus,                     ln_mus_ups,   &   ! MUSCL
+         &                 ln_trcadv_ubs,           nn_ubs_v,               &   ! UBS
+         &                 ln_trcadv_qck                                        ! QCK
       !!----------------------------------------------------------------------
       !
@@ -165 +189 @@
          WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtrc_adv : chose a advection scheme for tracers'
-         WRITE(numout,*) '      2nd order advection scheme     ln_trcadv_cen2   = ', ln_trcadv_cen2
-         WRITE(numout,*) '      TVD advection scheme           ln_trcadv_tvd    = ', ln_trcadv_tvd
-         WRITE(numout,*) '      MUSCL  advection scheme        ln_trcadv_muscl  = ', ln_trcadv_muscl
-         WRITE(numout,*) '      MUSCL2 advection scheme        ln_trcadv_muscl2 = ', ln_trcadv_muscl2
-         WRITE(numout,*) '      UBS    advection scheme        ln_trcadv_ubs    = ', ln_trcadv_ubs
-         WRITE(numout,*) '      QUICKEST advection scheme      ln_trcadv_qck    = ', ln_trcadv_qck
+         WRITE(numout,*) '      centered scheme                           ln_trcadv_cen = ', ln_trcadv_cen
+         WRITE(numout,*) '            horizontal 2nd/4th order               nn_cen_h   = ', nn_fct_h
+         WRITE(numout,*) '            vertical   2nd/4th order               nn_cen_v   = ', nn_fct_v
+         WRITE(numout,*) '      Flux Corrected Transport scheme           ln_trcadv_fct = ', ln_trcadv_fct
+         WRITE(numout,*) '            horizontal 2nd/4th order               nn_fct_h   = ', nn_fct_h
+         WRITE(numout,*) '            vertical   2nd/4th order               nn_fct_v   = ', nn_fct_v
+         WRITE(numout,*) '            2nd order + vertical sub-timestepping  nn_fct_zts = ', nn_fct_zts
+         WRITE(numout,*) '      MUSCL scheme                              ln_trcadv_mus = ', ln_trcadv_mus
+         WRITE(numout,*) '            + upstream scheme near river mouths    ln_mus_ups = ', ln_mus_ups
+         WRITE(numout,*) '      UBS scheme                                ln_trcadv_ubs = ', ln_trcadv_ubs
+         WRITE(numout,*) '            vertical   2nd/4th order               nn_ubs_v   = ', nn_ubs_v
+         WRITE(numout,*) '      QUICKEST scheme                           ln_trcadv_qck = ', ln_trcadv_qck
       ENDIF
       !
 
       ioptio = 0                      ! Parameter control
-      IF( ln_trcadv_cen2   )   ioptio = ioptio + 1
-      IF( ln_trcadv_tvd    )   ioptio = ioptio + 1
-      IF( ln_trcadv_muscl  )   ioptio = ioptio + 1
-      IF( ln_trcadv_muscl2 )   ioptio = ioptio + 1
-      IF( ln_trcadv_ubs    )   ioptio = ioptio + 1
-      IF( ln_trcadv_qck    )   ioptio = ioptio + 1
-      !
+      IF( ln_trcadv_cen                      )   nadv = np_CEN
+      IF( ln_trcadv_fct                      )   nadv = np_FCT
+      IF( ln_trcadv_fct .AND. nn_fct_zts > 0 )   nadv = np_FCT_zts
+      IF( ln_trcadv_mus                      )   nadv = np_MUS
+      IF( ln_trcadv_ubs                      )   nadv = np_UBS
+      IF( ln_trcadv_qck                      )   nadv = np_QCK
+      !
+      IF( ioptio == 0 ) THEN
+         nadv = np_NO_adv
+         CALL ctl_warn( 'trc_adv_init: You are running without tracer advection.' )
+      ENDIF
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtrc_adv' )
       !
+      IF( ln_trcadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 )   &
+                        .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 )   ) THEN
+        CALL ctl_stop( 'trc_adv_init: CEN scheme, choose 2nd or 4th order' )
+      ENDIF
+      IF( ln_trcadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 )   &
+                        .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 )   ) THEN
+        CALL ctl_stop( 'trc_adv_init: FCT scheme, choose 2nd or 4th order' )
+      ENDIF
+      IF( ln_trcadv_fct .AND. nn_fct_zts > 0 ) THEN
+         IF( nn_fct_h == 4 ) THEN
+            nn_fct_h = 2
+            CALL ctl_stop( 'trc_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' )
+         ENDIF
+         IF( lk_vvl ) THEN
+            CALL ctl_stop( 'trc_adv_init: vertical sub-timestepping not allow in non-linear free surface' )
+         ENDIF
+         IF( nn_fct_zts == 1 )   CALL ctl_warn( 'trc_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' )
+      ENDIF
+      IF( ln_trcadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN
+        CALL ctl_stop( 'trc_adv_init: UBS scheme, choose 2nd or 4th order' )
+      ENDIF
+      IF( ln_trcadv_ubs .AND. nn_ubs_v == 4 ) THEN
+         CALL ctl_warn( 'trc_adv_init: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' )
+      ENDIF
+      IF( ln_isfcav ) THEN                                                       ! ice-shelf cavities
+         IF(  ln_trcadv_cen .AND. nn_cen_v /= 4    .OR.   &                            ! NO 4th order with ISF
+            & ln_trcadv_fct .AND. nn_fct_v /= 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' )
+      ENDIF
+      !
       !                              ! Set nadv
-      IF( ln_trcadv_cen2   )   nadv =  1
-      IF( ln_trcadv_tvd    )   nadv =  2
-      IF( ln_trcadv_muscl  )   nadv =  3
-      IF( ln_trcadv_muscl2 )   nadv =  4
-      IF( ln_trcadv_ubs    )   nadv =  5
-      IF( ln_trcadv_qck    )   nadv =  6
+      IF( ln_trcadv_cen  )   nadv = np_CEN
+      IF( ln_trcadv_fct  )   nadv = np_FCT
+      IF( nn_fct_zts > 0 )   nadv = np_FCT_zts
+      IF( ln_trcadv_mus  )   nadv = np_MUS
+      IF( ln_trcadv_ubs  )   nadv = np_UBS
+      IF( ln_trcadv_qck  )   nadv = np_QCK
       !
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used'
-         IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used'
-         IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used'
-         IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    scheme is used'
-         IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
-         IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
+         IF( nadv == np_NO_adv  )   WRITE(numout,*) '         NO passive tracer advection'
+         IF( nadv == np_CEN     )   WRITE(numout,*) '         CEN      scheme is used. Horizontal order: ', nn_cen_h,   &
+            &                                                                        ' Vertical   order: ', nn_cen_v
+         IF( nadv == np_FCT     )   WRITE(numout,*) '         FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
+            &                                                                       ' Vertical   order: ', nn_fct_v
+         IF( nadv == np_FCT_zts )   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping'
+         IF( nadv == np_MUS     )   WRITE(numout,*) '         MUSCL    scheme is used'
+         IF( nadv == np_UBS     )   WRITE(numout,*) '         UBS      scheme is used'
+         IF( nadv == np_QCK     )   WRITE(numout,*) '         QUICKEST scheme is used'
       ENDIF
       !

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90

@@ -73 +73 @@
       IF( nn_timing == 1 )   CALL timing_start('trc_ldf')
       !
-      
       IF( l_trdtrc )  THEN
          CALL wrk_alloc( jpi,jpj,jpk,jptra,   ztrtrd )
          ztrtrd(:,:,:,:)  = tra(:,:,:,:)
       ENDIF
-
-      !                                        ! set the lateral diffusivity coef. for passive tracer      
+      !
+      !                                        !* set the lateral diffusivity coef. for passive tracer      
       CALL wrk_alloc( jpi,jpj,jpk,   zahu, zahv )
       zahu(:,:,:) = rldf * ahtu(:,:,:)
@@ -86 +85 @@
       SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend
       !
-      CASE ( n_lap   )                                ! iso-level laplacian
+      CASE ( np_lap   )                               ! iso-level laplacian
          CALL tra_ldf_lap  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb,      tra, jptra,  1   )
          !
-      CASE ( n_lap_i )                                ! laplacian : standard iso-neutral operator (Madec)
+      CASE ( np_lap_i )                               ! laplacian : standard iso-neutral operator (Madec)
          CALL tra_ldf_iso  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,  1   )
          !
-      CASE ( n_lap_it )                               ! laplacian : triad iso-neutral operator (griffies)
+      CASE ( np_lap_it )                              ! laplacian : triad iso-neutral operator (griffies)
          CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,  1   )
          !
-      CASE ( n_blp , n_blp_i , n_blp_it )             ! bilaplacian: all operator (iso-level, -neutral)
+      CASE ( np_blp , np_blp_i , np_blp_it )          ! bilaplacian: all operator (iso-level, -neutral)
          CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb     , tra, jptra, nldf )
          !
@@ -119 +118 @@
    END SUBROUTINE trc_ldf
 
-!!gm ldf_ctl should be called in trcini  so that l_ldfslp=T  cause the slope init and calculation
 
    SUBROUTINE trc_ldf_ini
@@ -128 +126 @@
       !!
       !! ** Method  :   set nldf from the namtra_ldf logicals
-      !!      nldf == -2   No lateral diffusion
       !!      nldf ==  0   laplacian operator
       !!      nldf ==  1   Rotated laplacian operator
@@ -134 +131 @@
       !!      nldf ==  3   Rotated bilaplacian
       !!----------------------------------------------------------------------
-      INTEGER ::   ioptio, ierr         ! temporary integers
-      INTEGER ::  ios                 ! Local integer output status for namelist read
-      !
-      NAMELIST/namtrc_ldf/ ln_trcldf_lap, ln_trcldf_blp, &
+      INTEGER ::   ioptio, ierr   ! temporary integers
+      INTEGER ::   ios            ! Local integer output status for namelist read
+      !
+      NAMELIST/namtrc_ldf/ ln_trcldf_lap, ln_trcldf_blp,                                  &
          &                 ln_trcldf_lev, ln_trcldf_hor, ln_trcldf_iso, ln_trcldf_triad,  &
          &                 rn_ahtrc_0   , rn_bhtrc_0
@@ -172 +169 @@
       IF( ln_trcldf_lap )   ioptio = ioptio + 1
       IF( ln_trcldf_blp )   ioptio = ioptio + 1
-      IF( ioptio >  1 )   CALL ctl_stop( 'trc_ldf_ctl: use ONE or NONE of the 2 lap/bilap operator type on tracer' )
-      IF( ioptio == 0 )   nldf = n_no_ldf   ! No lateral diffusion
+      IF( ioptio >  1   )   CALL ctl_stop( 'trc_ldf_ctl: use ONE or NONE of the 2 lap/bilap operator type on tracer' )
+      IF( ioptio == 0   )   nldf = np_no_ldf   ! No lateral diffusion
       
       IF( ln_trcldf_lap .AND. ln_trcldf_blp )   CALL ctl_stop( 'trc_ldf_ctl: bilaplacian should be used on both TRC and TRA' )
@@ -189 +186 @@
       IF( ln_trcldf_lap ) THEN      !==  laplacian operator  ==!
          IF ( ln_zco ) THEN                ! z-coordinate
-            IF ( ln_trcldf_lev   )   nldf = n_lap     ! iso-level = horizontal (no rotation)
-            IF ( ln_trcldf_hor   )   nldf = n_lap     ! iso-level = horizontal (no rotation)
-            IF ( ln_trcldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard  (   rotation)
-            IF ( ln_trcldf_triad )   nldf = n_lap_it  ! iso-neutral: triad     (   rotation)
+            IF ( ln_trcldf_lev   )   nldf = np_lap     ! iso-level = horizontal (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = np_lap     ! iso-level = horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation)
+            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation)
          ENDIF
          IF ( ln_zps ) THEN             ! z-coordinate with partial step
             IF ( ln_trcldf_lev   )   ierr = 1         ! iso-level not allowed 
-            IF ( ln_trcldf_hor   )   nldf = n_lap     ! horizontal (no rotation)
-            IF ( ln_trcldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation)
-            IF ( ln_trcldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_trcldf_hor   )   nldf = np_lap     ! horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard (rotation)
+            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad    (rotation)
          ENDIF
          IF ( ln_sco ) THEN             ! s-coordinate
-            IF ( ln_trcldf_lev   )   nldf = n_lap     ! iso-level  (no rotation)
-            IF ( ln_trcldf_hor   )   nldf = n_lap_it  ! horizontal (   rotation)       !!gm   a checker....
-            IF ( ln_trcldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation)
-            IF ( ln_trcldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_trcldf_lev   )   nldf = np_lap     ! iso-level  (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = np_lap_it  ! horizontal (   rotation)       !!gm   a checker....
+            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard (rotation)
+            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad    (rotation)
          ENDIF
          !                                ! diffusivity ratio: passive / active tracers 
@@ -217 +214 @@
          ENDIF
       ENDIF
-
+      !
       IF( ln_trcldf_blp ) THEN      !==  bilaplacian operator  ==!
          IF ( ln_zco ) THEN                ! z-coordinate
-            IF ( ln_trcldf_lev   )   nldf = n_blp     ! iso-level = horizontal (no rotation)
-            IF ( ln_trcldf_hor   )   nldf = n_blp     ! iso-level = horizontal (no rotation)
-            IF ( ln_trcldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
-            IF ( ln_trcldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_trcldf_lev   )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = np_blp     ! iso-level = horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation)
+            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation)
          ENDIF
          IF ( ln_zps ) THEN             ! z-coordinate with partial step
             IF ( ln_trcldf_lev   )   ierr = 1         ! iso-level not allowed 
-            IF ( ln_trcldf_hor   )   nldf = n_blp     ! horizontal (no rotation)
-            IF ( ln_trcldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
-            IF ( ln_trcldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_trcldf_hor   )   nldf = np_blp     ! horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation)
+            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation)
          ENDIF
          IF ( ln_sco ) THEN             ! s-coordinate
-            IF ( ln_trcldf_lev   )   nldf = n_blp     ! iso-level  (no rotation)
-            IF ( ln_trcldf_hor   )   nldf = n_blp_it  ! horizontal (   rotation)       !!gm   a checker....
-            IF ( ln_trcldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
-            IF ( ln_trcldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+            IF ( ln_trcldf_lev   )   nldf = np_blp     ! iso-level  (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = np_blp_it  ! horizontal (   rotation)       !!gm   a checker....
+            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation)
+            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation)
          ENDIF
          !                                ! diffusivity ratio: passive / active tracers 
@@ -248 +245 @@
          ENDIF
       ENDIF
-
+      !
       IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
       IF( ln_ldfeiv .AND. .NOT.ln_trcldf_iso )   &
@@ -256 +253 @@
          IF( .NOT.l_ldfslp )   CALL ctl_stop( '          the rotation of the diffusive tensor require l_ldfslp' )
       ENDIF
-
+      !
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nldf == n_no_ldf )   WRITE(numout,*) '          NO lateral diffusion'
-         IF( nldf == n_lap    )   WRITE(numout,*) '          laplacian iso-level operator'
-         IF( nldf == n_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)'
-         IF( nldf == n_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)'
-         IF( nldf == n_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator'
-         IF( nldf == n_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)'
-         IF( nldf == n_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)'
+         IF( nldf == np_no_ldf )   WRITE(numout,*) '          NO lateral diffusion'
+         IF( nldf == np_lap    )   WRITE(numout,*) '          laplacian iso-level operator'
+         IF( nldf == np_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)'
+         IF( nldf == np_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)'
+         IF( nldf == np_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator'
+         IF( nldf == np_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)'
+         IF( nldf == np_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)'
       ENDIF
       !

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90

@@ -11 +11 @@
    !!   'key_top'                                                TOP models
    !!----------------------------------------------------------------------
-   !!   trc_zdf     : update the tracer trend with the lateral diffusion
-   !!   trc_zdf_ini : initialization, namelist read, and parameters control
+   !!   trc_zdf      : update the tracer trend with the lateral diffusion
+   !!   trc_zdf_ini  : initialization, namelist read, and parameters control
    !!----------------------------------------------------------------------
-   USE oce_trc         ! ocean dynamics and active tracers
-   USE trc             ! ocean passive tracers variables
-   USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine)
-   USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp     routine)
-   USE trcldf   
-   USE trd_oce
-   USE trdtra
-   USE prtctl_trc      ! Print control
+   USE trc           ! ocean passive tracers variables
+   USE oce_trc       ! ocean dynamics and active tracers
+   USE trd_oce       ! trends: ocean variables
+   USE trazdf_exp    ! vertical diffusion: explicit (tra_zdf_exp     routine)
+   USE trazdf_imp    ! vertical diffusion: implicit (tra_zdf_imp     routine)
+   USE trcldf        ! passive tracers: lateral diffusion
+   USE trdtra        ! trends manager: tracers 
+   USE prtctl_trc    ! Print control
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   trc_zdf          ! called by step.F90 
-   PUBLIC   trc_zdf_alloc    ! called by nemogcm.F90 
-   PUBLIC   trc_zdf_ini    ! called by nemogcm.F90 
-   !                                        !!: ** Vertical diffusion
-   !                                        (nam_trczdf) **
+   PUBLIC   trc_zdf         ! called by step.F90 
+   PUBLIC   trc_zdf_ini     ! called by nemogcm.F90 
+   PUBLIC   trc_zdf_alloc   ! called by nemogcm.F90 
+   
+   !                                        !!** Vertical diffusion (nam_trczdf) **
    LOGICAL , PUBLIC ::   ln_trczdf_exp       !: explicit vertical diffusion scheme flag
    INTEGER , PUBLIC ::   nn_trczdf_exp       !: number of sub-time step (explicit time stepping)
@@ -44 +44 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! NEMO/TOP 3.7 , NEMO Consortium (2015)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -59 +59 @@
       !
    END FUNCTION trc_zdf_alloc
+
 
    SUBROUTINE trc_zdf( kt )
@@ -87 +88 @@
 
       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra ) 
-         WRITE(charout, FMT="('zdf1 ')") ;  CALL prt_ctl_trc_info(charout)
-                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dt,                trb, tra, jptra ) 
-         WRITE(charout, FMT="('zdf2 ')") ;  CALL prt_ctl_trc_info(charout)
-                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
       CASE ( 0 ) ;  CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
       CASE ( 1 ) ;  CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme          
-
       END SELECT
 
@@ -118 +111 @@
    END SUBROUTINE trc_zdf
 
+
    SUBROUTINE trc_zdf_ini
       !!----------------------------------------------------------------------
@@ -132 +126 @@
       !!----------------------------------------------------------------------
       INTEGER ::  ios                 ! Local integer output status for namelist read
-      !
+      !!
       NAMELIST/namtrc_zdf/ ln_trczdf_exp  , nn_trczdf_exp
       !!----------------------------------------------------------------------
-
-      !                                ! Vertical mixing
+      !
       REWIND( numnat_ref )             ! namtrc_zdf in reference namelist 
       READ  ( numnat_ref, namtrc_zdf, IOSTAT = ios, ERR = 905)
 905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in reference namelist', lwp )
-
+      !
       REWIND( numnat_cfg )             ! namtrc_zdf in configuration namelist 
       READ  ( numnat_cfg, namtrc_zdf, IOSTAT = ios, ERR = 906 )
 906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in configuration namelist', lwp )
       IF(lwm) WRITE ( numont, namtrc_zdf )
-
-      IF(lwp) THEN                     !   ! Control print
+      !
+      IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
          WRITE(numout,*) '   Namelist namtrc_zdf : set vertical diffusion  parameters'
@@ -153 +146 @@
       ENDIF
 
-      !  Define the vertical tracer physics scheme
-      IF( ln_trczdf_exp ) THEN  ;  nzdf = 0     ! explicit scheme
-      ELSE                      ;  nzdf = 1     ! implicit scheme
+      !                                ! Define the vertical tracer physics scheme
+      IF( ln_trczdf_exp ) THEN   ;   nzdf = 0     ! explicit scheme
+      ELSE                       ;   nzdf = 1     ! implicit scheme
       ENDIF
 
-      ! Force implicit schemes
-      IF( ln_trcldf_iso                               )   nzdf = 1      ! iso-neutral lateral physics
-      IF( ln_trcldf_hor .AND. ln_sco                  )   nzdf = 1      ! horizontal lateral physics in s-coordinate
+      !                                ! Force implicit schemes
+      IF( ln_trcldf_iso              )   nzdf = 1      ! iso-neutral lateral physics
+      IF( ln_trcldf_hor .AND. ln_sco )   nzdf = 1      ! horizontal lateral physics in s-coordinate
 #if defined key_zdftke || defined key_zdfgls 
-                                                          nzdf = 1      ! TKE or GLS physics       
+                                         nzdf = 1      ! TKE or GLS physics       
 #endif
       IF( ln_trczdf_exp .AND. nzdf == 1 )  & 
@@ -175 +168 @@
          IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
       ENDIF
-
+      !
    END SUBROUTINE trc_zdf_ini
+   
 #else
    !!----------------------------------------------------------------------

branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/trcini.F90

@@ -53 +53 @@
       !!                or read data or analytical formulation
       !!---------------------------------------------------------------------
-      !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('trc_init')
@@ -92 +91 @@
       !
    END SUBROUTINE trc_init
+
 
    SUBROUTINE trc_ini_ctl
@@ -103 +103 @@
       l_trcdm2dc = ln_dm2dc .OR. ( ln_cpl .AND. ncpl_qsr_freq /= 1 )
       l_trcdm2dc = l_trcdm2dc  .AND. .NOT. lk_offline
-      IF( l_trcdm2dc .AND. lwp ) &
-         &   CALL ctl_warn(' Coupling with passive tracers and used of diurnal cycle. &
-         & Computation of a daily mean shortwave for some biogeochemical models) ')
-
-      !
-      IF( nn_cla == 1 )   &
-         &  CALL ctl_stop( ' Cross Land Advection not yet implemented with passive tracer ; nn_cla must be 0' )
+      IF( l_trcdm2dc .AND. lwp )   CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.',   &
+         &                           'Computation of a daily mean shortwave for some biogeochemical models ' )
       !
    END SUBROUTINE trc_ini_ctl
+
 
    SUBROUTINE trc_ini_inv
@@ -157 +153 @@
    END SUBROUTINE trc_ini_inv
 
+
    SUBROUTINE trc_ini_sms
       !!----------------------------------------------------------------------
@@ -196 +193 @@
       !
    END SUBROUTINE trc_ini_trp
+
 
    SUBROUTINE trc_ini_state
@@ -248 +246 @@
    END SUBROUTINE trc_ini_state
 
+
    SUBROUTINE top_alloc
       !!----------------------------------------------------------------------