Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/eosbn2.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/eosbn2.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/eosbn2.F90	(revision 12143)
@@ -30,4 +30,6 @@
    !!   eos_insitu_2d : Compute the in situ density for 2d fields
    !!   bn2           : Compute the Brunt-Vaisala frequency
+   !!   bn2           : compute the Brunt-Vaisala frequency
+   !!   eos_pt_from_ct: compute the potential temperature from the Conservative Temperature
    !!   eos_rab       : generic interface of in situ thermal/haline expansion ratio 
    !!   eos_rab_3d    : compute in situ thermal/haline expansion ratio
@@ -74,7 +76,7 @@
 
    !                               !!** Namelist nameos **
-   LOGICAL , PUBLIC ::   ln_TEOS10   ! determine if eos_pt_from_ct is used to compute sst_m
-   LOGICAL , PUBLIC ::   ln_EOS80   ! determine if eos_pt_from_ct is used to compute sst_m
-   LOGICAL , PUBLIC ::   ln_SEOS   ! determine if eos_pt_from_ct is used to compute sst_m
+   LOGICAL , PUBLIC ::   ln_TEOS10
+   LOGICAL , PUBLIC ::   ln_EOS80
+   LOGICAL , PUBLIC ::   ln_SEOS
 
    ! Parameters
@@ -1235,9 +1237,9 @@
       REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state
       READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' )
       !
       REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state
       READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' )
       IF(lwm) WRITE( numond, nameos )
       !
@@ -1647,4 +1649,7 @@
          !
       CASE( np_seos )                        !==  Simplified EOS     ==!
+
+         r1_S0  = 0.875_wp/35.16504_wp   ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct)
+         
          IF(lwp) THEN
             WRITE(numout,*)
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traadv.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traadv.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traadv.F90	(revision 12143)
@@ -196,9 +196,9 @@
       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 )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_adv in reference namelist' )
       !
       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 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_adv in configuration namelist' )
       IF(lwm) WRITE( numond, namtra_adv )
       !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traadv_fct.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traadv_fct.F90	(revision 12143)
@@ -21,4 +21,5 @@
    USE diaar5         ! AR5 diagnostics
    USE phycst  , ONLY : rau0_rcp
+   USE zdf_oce , ONLY : ln_zad_Aimp
    !
    USE in_out_manager ! I/O manager
@@ -86,4 +87,6 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup
+      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection
       !!----------------------------------------------------------------------
       !
@@ -97,4 +100,5 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
+      ll_zAimp = .FALSE.
       IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
       IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
@@ -116,4 +120,22 @@
       !
       zwi(:,:,:) = 0._wp        
+      !
+      ! If adaptive vertical advection, check if it is needed on this PE at this time
+      IF( ln_zad_Aimp ) THEN
+         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.
+      END IF
+      ! If active adaptive vertical advection, build tridiagonal matrix
+      IF( ll_zAimp ) THEN
+         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
+                  zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk)
+                  zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t_a(ji,jj,jk)
+                  zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      END IF
       !
       DO jn = 1, kjpt            !==  loop over the tracers  ==!
@@ -169,4 +191,29 @@
             END DO
          END DO
+         
+         IF ( ll_zAimp ) THEN
+            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
+            !
+            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     ztw(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes
+                  END DO
+               END DO
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
+                        &                                  * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+         END IF
          !                
          IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
@@ -277,4 +324,30 @@
             zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
          ENDIF
+         !         
+         IF ( ll_zAimp ) THEN
+            DO jk = 1, jpkm1     !* trend and after field with monotonic scheme
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     !                             ! total intermediate advective trends
+                     ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                        &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                        &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+                     ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            !
+            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
+            !
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF
          !
          CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. )
@@ -289,11 +362,35 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.  
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
-                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
-                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) &
-                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
+                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
+                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk)
+                  zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF ( ll_zAimp ) THEN
+            !
+            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
+            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
+                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
+                     ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic
+                  END DO
+               END DO
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.  
+                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) &
+                        &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF         
          !
          IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport
@@ -318,4 +415,7 @@
       END DO                     ! end of tracer loop
       !
+      IF ( ll_zAimp ) THEN
+         DEALLOCATE( zwdia, zwinf, zwsup )
+      ENDIF
       IF( l_trd .OR. l_hst ) THEN 
          DEALLOCATE( ztrdx, ztrdy, ztrdz )
@@ -559,5 +659,5 @@
          DO ji = fs_2, fs_jpim1
             ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point
-            ikb = mbkt(ji,jj)                !     -   above the last wet point
+            ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point
             !
             zwd (ji,jj,ikt) = 1._wp          ! top
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbc.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbc.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbc.F90	(revision 12143)
@@ -135,9 +135,9 @@
       REWIND( numnam_ref )              ! Namelist nambbc in reference namelist : Bottom momentum boundary condition
       READ  ( numnam_ref, nambbc, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbc in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbc in reference namelist' )
       !
       REWIND( numnam_cfg )              ! Namelist nambbc in configuration namelist : Bottom momentum boundary condition
       READ  ( numnam_cfg, nambbc, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbc in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbc in configuration namelist' )
       IF(lwm) WRITE ( numond, nambbc )
       !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbl.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbl.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trabbl.F90	(revision 12143)
@@ -485,9 +485,9 @@
       REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
       READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist' )
       !
       REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
       READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist' )
       IF(lwm) WRITE ( numond, nambbl )
       !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tradmp.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tradmp.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tradmp.F90	(revision 12143)
@@ -179,9 +179,9 @@
       REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation
       READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' )
       !
       REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation
       READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' )
       IF(lwm) WRITE ( numond, namtra_dmp )
       !
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traldf_iso.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traldf_iso.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traldf_iso.F90	(revision 12143)
@@ -289,5 +289,5 @@
          !!----------------------------------------------------------------------
          !
-         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp
+         ztfw(fs_2:1,:,:) = 0._wp     ;     ztfw(jpi:fs_jpim1,:,:) = 0._wp   ! avoid to potentially manipulate NaN values
          !
          ! Vertical fluxes
@@ -323,5 +323,5 @@
          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
             DO jk = 2, jpkm1       
-               DO jj = 1, jpjm1
+               DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   &
@@ -336,5 +336,5 @@
             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
                DO jk = 2, jpkm1 
-                  DO jj = 1, jpjm1
+                  DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
@@ -346,5 +346,5 @@
             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
                DO jk = 2, jpkm1 
-                  DO jj = 1, jpjm1
+                  DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      &
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tramle.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tramle.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tramle.F90	(revision 12143)
@@ -268,9 +268,9 @@
       REWIND( numnam_ref )              ! Namelist namtra_mle in reference namelist : Tracer advection scheme
       READ  ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_mle in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_mle in reference namelist' )
 
       REWIND( numnam_cfg )              ! Namelist namtra_mle in configuration namelist : Tracer advection scheme
       READ  ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_mle in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_mle in configuration namelist' )
       IF(lwm) WRITE ( numond, namtra_mle )
 
Index: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traqsr.F90	(revision 12062)
+++ NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/traqsr.F90	(revision 12143)
@@ -168,5 +168,5 @@
                DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
                   DO ji = fs_2, fs_jpim1
-                     zchl    = sf_chl(1)%fnow(ji,jj,1)
+                     zchl    = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
                      zCtot   = 40.6  * zchl**0.459
                      zze     = 568.2 * zCtot**(-0.746)
@@ -338,9 +338,9 @@
       REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist
       READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
       !
       REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist
       READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
       IF(lwm) WRITE ( numond, namtra_qsr )
       !
