Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/abl.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/abl.F90	(revision 12587)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/abl.F90	(revision 12588)
@@ -29,5 +29,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avm_abl      !: turbulent viscosity   [m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avt_abl      !: turbulent diffusivity [m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxl_abl      !: mixing length         [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxld_abl     !: dissipative mixing length    [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxlm_abl     !: master mixing length         [m]
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   tke_abl      !: turbulent kinetic energy [m2/s2]
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   fft_abl      !: Coriolis parameter    [1/s]
@@ -38,4 +39,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   msk_abl
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   rest_eq
+   
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::    cft_abl
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:  )     ::   taux_abl
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:  )     ::   tauy_abl  
    !
    INTEGER , PUBLIC :: nt_n, nt_a       !: now / after indices (equal 1 or 2)
@@ -55,16 +60,21 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( u_abl  (1:jpi,1:jpj,1:jpka,jptime), &
-         &      v_abl  (1:jpi,1:jpj,1:jpka,jptime), &
-         &      tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), &
-         &      avm_abl(1:jpi,1:jpj,1:jpka), &
-         &      avt_abl(1:jpi,1:jpj,1:jpka), &
-         &      mxl_abl(1:jpi,1:jpj,1:jpka), &         
-         &      tke_abl(1:jpi,1:jpj,1:jpka,jptime), &
-         &      fft_abl(1:jpi,1:jpj), &
-         &      pblh   (1:jpi,1:jpj), &         
-         &      msk_abl(1:jpi,1:jpj), &
-         &      rest_eq(1:jpi,1:jpj), &         
-         &      e3t_abl(1:jpka), e3w_abl(1:jpka), ght_abl(1:jpka), ghw_abl(1:jpka),                 STAT=ierr )
+      ALLOCATE( u_abl   (1:jpi,1:jpj,1:jpka,jptime     ), &
+         &      v_abl   (1:jpi,1:jpj,1:jpka,jptime     ), &
+         &      tq_abl  (1:jpi,1:jpj,1:jpka,jptime,jptq), &
+         &      tke_abl (1:jpi,1:jpj,1:jpka,jptime     ), &
+         &      avm_abl (1:jpi,1:jpj,1:jpka            ), &
+         &      avt_abl (1:jpi,1:jpj,1:jpka            ), &
+         &      mxld_abl(1:jpi,1:jpj,1:jpka            ), &         
+         &      mxlm_abl(1:jpi,1:jpj,1:jpka            ), & 
+         &      cft_abl (1:jpi,1:jpj,1:jpka            ), & 
+         &      fft_abl (1:jpi,1:jpj                   ), &
+         &      pblh    (1:jpi,1:jpj                   ), &         
+         &      taux_abl(1:jpi,1:jpj                   ), &   
+         &      tauy_abl(1:jpi,1:jpj                   ), &   		 
+         &      msk_abl (1:jpi,1:jpj                   ), &
+         &      rest_eq (1:jpi,1:jpj                   ), &         
+         &      e3t_abl (1:jpka), e3w_abl(1:jpka)       , &
+         &      ght_abl (1:jpka), ghw_abl(1:jpka)       , STAT=ierr )
          !
       abl_alloc = ierr
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90	(revision 12587)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90	(revision 12588)
@@ -2,10 +2,10 @@
    !!======================================================================
    !!                       ***  MODULE  ablmod  ***
-   !! Surface module :  ABL computation to provide atmospheric data 
+   !! Surface module :  ABL computation to provide atmospheric data
    !!                   for surface fluxes computation
    !!======================================================================
    !! History :  3.6  ! 2019-03  (F. Lemarié & G. Samson)  Original code
    !!----------------------------------------------------------------------
-   
+
    !!----------------------------------------------------------------------
    !!   abl_stp       : ABL single column model
@@ -16,5 +16,5 @@
 
    USE phycst         ! physical constants
-   USE dom_oce, ONLY  : tmask  
+   USE dom_oce, ONLY  : tmask
    USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa
    USE sbcblk         ! use rn_?fac
@@ -30,5 +30,5 @@
 
    PUBLIC   abl_stp   ! called by sbcabl.F90
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ustar2
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ustar2, zrough
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -39,14 +39,15 @@
 !===================================================================================================
    SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in
-              &            pu_dta, pv_dta, pt_dta, pq_dta,    & 
+              &            pu_dta, pv_dta, pt_dta, pq_dta,    &
               &            pslp_dta, pgu_dta, pgv_dta,        &
               &            pcd_du, psen, pevp,                &                            ! in/out
-              &            pwndm, ptaui, ptauj, ptaum         & 
-#if defined key_si3			  
-              &          , ptm_su,pssu_ice,pssv_ice,pssq_ice,pcd_du_ice  & 
-              &          , psen_ice, pevp_ice, pwndm_ice, pfrac_oce      & 
-              &          , ptaui_ice, ptauj_ice                          &
-#endif			  
-			  &   )
+              &            pwndm, ptaui, ptauj, ptaum         &
+#if defined key_si3
+              &          , ptm_su, pssu_ice, pssv_ice         &
+              &          , pssq_ice, pcd_du_ice, psen_ice     &
+              &          , pevp_ice, pwndm_ice, pfrac_oce     &
+              &          , ptaui_ice, ptauj_ice               &
+#endif
+              &      )
 !---------------------------------------------------------------------------------------------------
 
@@ -54,9 +55,9 @@
       !!                    ***  ROUTINE abl_stp ***
       !!
-      !! ** Purpose :   Time-integration of the ABL model 
+      !! ** Purpose :   Time-integration of the ABL model
       !!
-      !! ** Method  :   Compute atmospheric variables : vertical turbulence 
+      !! ** Method  :   Compute atmospheric variables : vertical turbulence
       !!                             + Coriolis term + newtonian relaxation
-      !!                
+      !!
       !! ** Action  : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
       !!              - Advance tracers to time n+1 (Euler backward scheme)
@@ -70,5 +71,5 @@
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psst       ! sea-surface temperature [Celsius]
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu       ! sea-surface u (U-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point) 
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point)
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq       ! sea-surface humidity
       REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pu_dta     ! large-scale windi
@@ -82,27 +83,25 @@
       REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   psen       ! Ch x Du
       REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pevp       ! Ce x Du
-      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd||      
+      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd||
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui      ! taux
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj      ! tauy
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau|| 
-      !
-#if defined key_si3	  
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau||
+      !
+#if defined key_si3
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptm_su       ! ice-surface temperature [K]
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu_ice     ! ice-surface u (U-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point) 	  
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity 
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du_ice   ! Cd x Du over ice (T-point)
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psen_ice     ! Ch x Du over ice (T-point)
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||  
-      !REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pfrac_oce	   !!GS: out useless ?
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce	   !
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce    ! ocean fraction
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (U-point)
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point) 	  
-#endif	  
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point)
+#endif
       !
       REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zwnd_i, zwnd_j
-      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF    
-      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)    ::   z_cft      !--FL--to be removed after the test phase   
+      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF
       !
       REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a
@@ -113,6 +112,7 @@
       REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
       REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
-      !
-      !!---------------------------------------------------------------------      
+      LOGICAL             ::   SemiImp_Cor = .FALSE.
+      !
+      !!---------------------------------------------------------------------
       !
       IF(lwp .AND. kt == nit000) THEN                  ! control print
@@ -120,18 +120,20 @@
          WRITE(numout,*) 'abl_stp : ABL time stepping'
          WRITE(numout,*) '~~~~~~'
-      ENDIF 
+      ENDIF
       !
       IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
-      !! Compute ustar squared as Cd || Uatm-Uoce ||^2 
-      !! needed for surface boundary condition of TKE 
+      IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) )
+      !! Compute ustar squared as Cd || Uatm-Uoce ||^2
+      !! needed for surface boundary condition of TKE
       !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
       DO_2D_11_11
          zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj)
 #if defined key_si3
-         zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 
-         ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 
+         zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)
+         ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice
 #else
-         ustar2(ji,jj) = zzoce	
+         ustar2(ji,jj) = zzoce
 #endif
+         zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT(Cdn_oce(ji,jj)) ) !<-- recover the value of z0 from Cdn_oce
       END_2D
       !
@@ -140,36 +142,37 @@
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj) 
-         
+      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
+
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  2 *** Advance tracers to time n+1
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      
+
       !-------------
       DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj)
-      !-------------    
-         ! Compute matrix elements for interior points 
+      !-------------
+         ! Compute matrix elements for interior points
          DO jk = 3, jpkam1
             DO ji = 1, jpi   ! vector opt.
-               z_elem_a( ji,     jk              ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
-               z_elem_c( ji,     jk              ) = - rDt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal       
-               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal           
-            END DO
-         END DO  
-         ! Boundary conditions  
-         DO ji = 1, jpi   ! vector opt. 
-            ! Neumann at the bottom          
-            z_elem_a( ji,     2              ) = 0._wp
-            z_elem_c( ji,     2              ) = - rDt_abl * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                         
+               z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
+               z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal
+               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal
+            END DO
+         END DO
+         ! Boundary conditions
+         DO ji = 1, jpi   ! vector opt.
+            ! Neumann at the bottom
+            z_elem_a( ji, 2 ) = 0._wp
+            z_elem_c( ji, 2    ) = - rDt_abl * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )
             ! Homogeneous Neumann at the top
-            z_elem_a( ji,     jpka           ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
-            z_elem_c( ji,     jpka           ) = 0._wp
-            z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
-         END DO            
+            z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )
+            z_elem_c( ji, jpka ) = 0._wp
+            z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
+         END DO
 
          DO jtra = 1,jptq  ! loop on active tracers
-               
+
             DO jk = 3, jpkam1
-               DO ji = 1,jpi
+               DO ji = 2, jpim1
+               !DO ji = 1,jpi  !!GS: to be checked
                   tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
                END DO
@@ -177,39 +180,39 @@
 
             IF(jtra == jp_ta) THEN
-               DO ji = 1,jpi  ! boundary conditions for temperature              
-                  zztmp1 = psen(ji, jj) 
-                  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )	
-#if defined key_si3				  
-                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 
+               DO ji = 1,jpi  ! surface boundary condition for temperature
+                  zztmp1 = psen(ji, jj)
+                  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )
+#if defined key_si3
+                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)
                   zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
-#endif				  
-                  z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1				  
-                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rDt_abl * zztmp2               
+#endif
+                  z_elem_b( ji,     2                ) = e3t_abl( 2    ) - z_elem_c( ji, 2                    ) + rDt_abl * zztmp1
+                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rDt_abl * zztmp2
                   tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
-               END DO 
-            ELSE   
-               DO ji = 1,jpi  ! boundary conditions for humidity              
-                  zztmp1 = pevp(ji, jj) 
-                  zztmp2 = pevp(ji, jj) * pssq(ji, jj)   
-#if defined key_si3				  
-                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 
-                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)    
-#endif	
-                  z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
-                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rDt_abl * zztmp2               
+               END DO
+            ELSE ! jp_qa
+               DO ji = 1,jpi  ! surface boundary condition for humidity
+                  zztmp1 = pevp(ji, jj)
+                  zztmp2 = pevp(ji, jj) * pssq(ji, jj)
+#if defined key_si3
+                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj)
+                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)
+#endif
+                  z_elem_b( ji,     2                ) = e3t_abl( 2    ) - z_elem_c( ji, 2                    ) + rDt_abl * zztmp1
+                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj,    2, nt_n, jtra ) + rDt_abl * zztmp2
                   tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
-               END DO                                     
+               END DO
             END IF
             !!
             !! Matrix inversion
             !! ----------------------------------------------------------
-            DO ji = 1,jpi            
+            DO ji = 1,jpi
                zcff                       =  1._wp / z_elem_b( ji, 2 )
                zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 )
                tq_abl(ji,jj,2,nt_a,jtra) =    zcff * tq_abl(ji,jj,2,nt_a,jtra)
-            END DO             
-
-            DO jk = 3, jpka            
-               DO ji = 1,jpi            
+            END DO
+
+            DO jk = 3, jpka
+               DO ji = 1,jpi
                   zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
                   zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
@@ -218,123 +221,168 @@
                END DO
             END DO
-!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...         
-            DO jk = jpkam1,2,-1            
-               DO ji = 1,jpi  
+            !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...
+            DO jk = jpkam1,2,-1
+               DO ji = 1,jpi
                   tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    &
                      &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
                END DO
             END DO
-         
-         END DO   !<-- loop on tracers                    
-         !! 
-      !-------------
-      END DO             ! end outer loop 
-      !-------------   
-      
-      
+
+         END DO   !<-- loop on tracers
+         !!
+      !-------------
+      END DO             ! end outer loop
+      !-------------
+
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  3 *** Compute Coriolis term with geostrophic guide
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
-      !-------------        
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( SemiImp_Cor ) THEN
+
+      !-------------
       DO jk = 2, jpka    ! outer loop
       !-------------
-         !             
+         !
          ! Advance u_abl & v_abl to time n+1
          DO_2D_11_11
             zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2
-   
+
             u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  &
                &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    &
                &                 +  rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  &
                &                               / (1._wp + gamma_Cor*gamma_Cor*zcff)
-               
+
             v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  &
                &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   &
                &                 -  rDt_abl * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) &
-               &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)                
+               &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)
          END_2D
-         !                                   
+         !
       !-------------
       END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
-      !-------------                
-      !
-      IF( ln_geos_winds ) THEN
-         DO jj = 1, jpj    ! outer loop       
+      !-------------
+      !
+      !IF( ln_geos_winds ) THEN
+      !   DO jj = 1, jpj    ! outer loop
+      !      DO jk = 1, jpka
+      !         DO ji = 1, jpi
+      !            u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
+      !               &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
+      !            v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
+      !               &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
+      !         END DO
+      !      END DO
+      !   END DO
+      !END IF
+
+      ELSE
+
+      !-------------
+      DO jk = 2, jpka    ! outer loop
+      !-------------
+         !
+         IF( MOD( kt, 2 ) == 0 ) then
+                 ! Advance u_abl & v_abl to time n+1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_n ) - pgv_dta(ji  ,jj  ,jk)  )
+               u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff
+               zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_a ) - pgu_dta(ji  ,jj  ,jk)  )
+               v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff )
+               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  u_abl( ji, jj, jk, nt_a )
+            END DO
+         END DO
+         !
+         ELSE
+         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_n ) - pgu_dta(ji  ,jj  ,jk)  )
+               v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff
+               zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_a ) - pgv_dta(ji  ,jj  ,jk)  )
+               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff )
+               v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  v_abl( ji, jj, jk, nt_a )
+            END DO
+         END DO
+         !
+         ENDIF
+         !
+      !-------------
+      END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
+      !-------------
+
+      ENDIF
+
+      !-------------
+      !
+      IF( ln_hpgls_frc ) THEN
+         DO jj = 1, jpj    ! outer loop
             DO jk = 1, jpka
-               DO ji = 1, jpi 
-                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
-                     &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
-                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
-                     &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
-               END DO
-            END DO
-         END DO      
-      END IF 
-      !-------------            
-      !
-      IF( ln_hpgls_frc ) THEN
-         DO jj = 1, jpj    ! outer loop       
-            DO jk = 1, jpka
-               DO ji = 1, jpi 
+               DO ji = 1, jpi
                   u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
-                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)  
+                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)
                ENDDO
             ENDDO
-         ENDDO 
+         ENDDO
       END IF
-	  
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                            !  4 *** Advance u,v to time n+1 
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
-      !
-      !  Vertical diffusion for u_abl     
+
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  4 *** Advance u,v to time n+1
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      !  Vertical diffusion for u_abl
       !-------------
       DO jj = 1, jpj    ! outer loop
-      !-------------    
+      !-------------
 
          DO jk = 3, jpkam1
-            DO ji = 1, jpi  
+            DO ji = 1, jpi
                z_elem_a( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal
-               z_elem_c( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal                
+               z_elem_c( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal
                z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal
             END DO
-         END DO  
-            
-         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)            
+         END DO
+
+         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)
             !++ Surface boundary condition
             z_elem_a( ji,     2    ) = 0._wp
-            z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
-            !
-			zztmp1  = pcd_du(ji, jj)
-			zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )			
-#if defined key_si3				  
+            z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )
+            !
+            zztmp1  = pcd_du(ji, jj)
+            zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) * rn_vfac
+#if defined key_si3
             zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
-			zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) )
-			zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
-#endif           
-			z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 			
+            zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) * rn_vfac
+            zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
+#endif
+            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
             u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rDt_abl * zztmp2
-			
-            !++ Top Neumann B.C.
-            !z_elem_a( ji,     jpka ) = - 0.5_wp * rDt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka ) 
-            !z_elem_c( ji,     jpka ) = 0._wp
-            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
-            !++ Top Dirichlet B.C.
-            z_elem_a( ji,     jpka )       = 0._wp
-            z_elem_c( ji,     jpka )       = 0._wp
-            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
-            u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)              
-         END DO 
+
+            IF( ln_topbc_neumann ) THEN
+               !++ Top Neumann B.C.
+               z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
+               z_elem_c( ji,     jpka ) = 0._wp
+               z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
+               !u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl   ( ji, jj, jpka, nt_a )
+            ELSE
+               !++ Top Dirichlet B.C.
+               z_elem_a( ji,     jpka )       = 0._wp
+               z_elem_c( ji,     jpka )       = 0._wp
+               z_elem_b( ji,     jpka )       = e3t_abl( jpka )
+               u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)
+            ENDIF
+
+         END DO
          !!
          !! Matrix inversion
          !! ----------------------------------------------------------
-         DO ji = 2, jpi          
+         DO ji = 2, jpi
             zcff                 =   1._wp / z_elem_b( ji, 2 )
-            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 ) 
+            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 )
             u_abl (ji,jj,2,nt_a) =    zcff * u_abl(ji,jj,2,nt_a)
-         END DO             
-
-         DO jk = 3, jpka            
-            DO ji = 2, jpi              
+         END DO
+
+         DO jk = 3, jpka
+            DO ji = 2, jpi
                zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
                zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
@@ -343,64 +391,70 @@
             END DO
          END DO
-            
-         DO jk = jpkam1,2,-1            
+
+         DO jk = jpkam1,2,-1
             DO ji = 2, jpi
                u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a)
             END DO
          END DO
-                                          
-      !-------------
-      END DO             ! end outer loop 
-      !-------------   
-
-      !
-      !  Vertical diffusion for v_abl     
+
+      !-------------
+      END DO             ! end outer loop
+      !-------------
+
+      !
+      !  Vertical diffusion for v_abl
       !-------------
       DO jj = 2, jpj   ! outer loop
-      !-------------    
+      !-------------
          !
          DO jk = 3, jpkam1
-            DO ji = 1, jpi   
+            DO ji = 1, jpi
                z_elem_a( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
-               z_elem_c( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal              
+               z_elem_c( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal
                z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal
             END DO
          END DO
 
-         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)          
+         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)
             !++ Surface boundary condition
             z_elem_a( ji,     2    ) = 0._wp
-            z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )        
-            !
-			zztmp1 = pcd_du(ji, jj)
-			zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) )  
-#if defined key_si3				  
+            z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )
+            !
+            zztmp1 = pcd_du(ji, jj)
+            zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) * rn_vfac
+#if defined key_si3
             zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
-			zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) )
-			zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
-#endif 			
-            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 
+            zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) ) * rn_vfac
+            zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
+#endif
+            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
             v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
-            !++ Top Neumann B.C.            
-            !z_elem_a( ji,     jpka ) = -rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
-            !z_elem_c( ji,     jpka ) = 0._wp
-            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
-            !++ Top Dirichlet B.C.
-            z_elem_a( ji,     jpka )       = 0._wp
-            z_elem_c( ji,     jpka )       = 0._wp
-            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
-            v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)             
-         END DO 
+
+            IF( ln_topbc_neumann ) THEN
+               !++ Top Neumann B.C.
+               z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )
+               z_elem_c( ji,     jpka ) = 0._wp
+               z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
+               !v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl   ( ji, jj, jpka, nt_a )
+            ELSE
+               !++ Top Dirichlet B.C.
+               z_elem_a( ji,     jpka )       = 0._wp
+               z_elem_c( ji,     jpka )       = 0._wp
+               z_elem_b( ji,     jpka )       = e3t_abl( jpka )
+               v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)
+            ENDIF
+
+         END DO
          !!
          !! Matrix inversion
          !! ----------------------------------------------------------
-         DO ji = 1, jpi              
+         DO ji = 1, jpi
             zcff                 =  1._wp / z_elem_b( ji, 2 )
-            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       ) 
-            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a ) 
-         END DO             
-
-         DO jk = 3, jpka            
-            DO ji = 1, jpi                   
+            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       )
+            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a )
+         END DO
+
+         DO jk = 3, jpka
+            DO ji = 1, jpi
                zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
                zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
@@ -409,37 +463,35 @@
             END DO
          END DO
-            
-         DO jk = jpkam1,2,-1            
-            DO ji = 1, jpi    
+
+         DO jk = jpkam1,2,-1
+            DO ji = 1, jpi
                v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a)
             END DO
          END DO
-         !                                          
-      !-------------
-      END DO             ! end outer loop 
-      !-------------   
-    
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                            !  5 *** Apply nudging on the dynamics and the tracers 
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
-      z_cft(:,:,:) = 0._wp      
-      
+         !
+      !-------------
+      END DO             ! end outer loop
+      !-------------
+
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  5 *** Apply nudging on the dynamics and the tracers
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
       IF( nn_dyn_restore > 0  ) THEN
-         !------------- 
+         !-------------
          DO jk = 2, jpka    ! outer loop
-         !-------------       
+         !-------------
             DO_2D_01_01
                zcff1 = pblh( ji, jj )
-               zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                        
-               zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
+               zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
+               zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )
                zmsk  = msk_abl(ji,jj)
                zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   &
                   &  + jp_alp1_dyn * zsig    + jp_alp0_dyn
                zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
-                                                             ! rn_Dt = rDt_abl / nn_fsbc                          
+                                                             ! rn_Dt = rDt_abl / nn_fsbc
                zcff  = zcff * rest_eq(ji,jj)
-               z_cft( ji, jj, jk ) = zcff
                u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   &
-                  &                               + zcff   * pu_dta( ji, jj, jk       )                      
+                  &                               + zcff   * pu_dta( ji, jj, jk       )
                v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   &
                   &                               + zcff   * pv_dta( ji, jj, jk       )
@@ -447,37 +499,37 @@
          !-------------
          END DO             ! end outer loop
-         !-------------               
+         !-------------
       END IF
 
-      !------------- 
+      !-------------
       DO jk = 2, jpka    ! outer loop
-      !-------------       
+      !-------------
          DO_2D_11_11
             zcff1 = pblh( ji, jj )
             zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
-            zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
+            zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )
             zmsk  = msk_abl(ji,jj)
             zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   &
                &  + jp_alp1_tra * zsig    + jp_alp0_tra
             zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points
-                                                          ! rn_Dt = rDt_abl / nn_fsbc                          
-            !z_cft( ji, jj, jk ) = zcff
+                                                          ! rn_Dt = rDt_abl / nn_fsbc
             tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   &
                &                                       + zcff   * pt_dta( ji, jj, jk )
-            
+
             tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   &
                &                                       + zcff   * pq_dta( ji, jj, jk )
-            
+
          END_2D
       !-------------
       END DO             ! end outer loop
-      !------------- 	  
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                            !  6 *** MPI exchanges  
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !
-      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a      ), 'T', -1. )
+      !-------------
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  6 *** MPI exchanges
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a)      , 'T', -1., kfillmode = jpfillnothing )
       CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T',  1., tq_abl(:,:,:,nt_a,jp_qa), 'T',  1., kfillmode = jpfillnothing )   ! ++++ this should not be needed...
       !
+#if defined key_iomput
       ! first ABL level
       IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) )
@@ -497,5 +549,5 @@
       IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka           ) )
       IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka           ) )
-      IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxl_abl(:,:,2:jpka           ) )
+      IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxlm_abl(:,:,2:jpka          ) )
       IF ( iom_use("pblh"   ) ) CALL iom_put (  "pblh"  ,    pblh(:,:                  ) )
       ! debug (to be removed)
@@ -504,5 +556,4 @@
       IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta",  pt_dta(:,:,2:jpka) )
       IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta",  pq_dta(:,:,2:jpka) )
-      IF ( iom_use("coeft") ) CALL iom_put ( "coeft",   z_cft(:,:,2:jpka) )
       IF( ln_geos_winds ) THEN
          IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2           ) )
@@ -513,37 +564,36 @@
          IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)   )
       END IF
-      !
+#endif
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  7 *** Finalize flux computation
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
-
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
       DO_2D_11_11
-         ztemp             = tq_abl  ( ji, jj, 2, nt_a, jp_ta ) 
-         zhumi             = tq_abl  ( ji, jj, 2, nt_a, jp_qa ) 
-         !zcff              = pslp_dta( ji, jj ) /   &              !<-- At this point ztemp and zhumi should not be zero ...
-         !   &                        (  R_dry*ztemp * ( 1._wp + rctv0*zhumi )  )
-         zcff              = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
-         psen ( ji, jj )   =      cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )
-         pevp ( ji, jj )   = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) )
-         rhoa( ji, jj )   = zcff              
+         ztemp          = tq_abl  ( ji, jj, 2, nt_a, jp_ta )
+         zhumi          = tq_abl  ( ji, jj, 2, nt_a, jp_qa )
+         zcff           = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
+         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )   !GS: negative sign to respect aerobulk convention
+         pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) )
+         rhoa( ji, jj ) = zcff
       END_2D
-      
+
       DO_2D_01_01
-         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) )  
-         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) ) 
+         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji,jj) + pssu(ji-1,jj  ) )
+         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj) + pssv(ji  ,jj-1) )
       END_2D
-      ! 
+      !
       CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
       !
       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
       DO_2D_11_11
-         zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
-            &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj)  )  ! * msk_abl(ji,jj)
-         zztmp         = rhoa(ji,jj) * pcd_du(ji,jj)
-         
-         pwndm (ji,jj) =         zcff
-         ptaum (ji,jj) = zztmp * zcff
-         zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
-         zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+         zcff            = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
+            &                   + zwnd_j(ji,jj) * zwnd_j(ji,jj)   )  ! * msk_abl(ji,jj)
+         zztmp           = rhoa(ji,jj) * pcd_du(ji,jj)
+         pwndm (ji,jj)   =         zcff
+         ptaum (ji,jj)   = zztmp * zcff
+         zwnd_i(ji,jj)   = zztmp * zwnd_i(ji,jj)
+         zwnd_j(ji,jj)   = zztmp * zwnd_j(ji,jj)
+         taux_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_i(ji,jj)
+         tauy_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_j(ji,jj)
       END_2D
       ! ... utau, vtau at U- and V_points, resp.
@@ -574,8 +624,8 @@
          ! ------------------------------------------------------------ !
          DO_2D_00_00
-            
-            zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) )
-            zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
-   
+
+            zztmp1 = 0.5_wp * ( u_abl(ji+1,jj  ,2,nt_a) + u_abl(ji,jj,2,nt_a) )
+            zztmp2 = 0.5_wp * ( v_abl(ji  ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
+
             ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
                &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          &
@@ -592,8 +642,8 @@
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  8 *** Swap time indices for the next timestep
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       nt_n = 1 + MOD( kt  , 2)
       nt_a = 1 + MOD( kt+1, 2)
-      !	  
+      !
 !---------------------------------------------------------------------------------------------------
    END SUBROUTINE abl_stp
@@ -634,13 +684,13 @@
       !!                (= Kz dz[Ub] * dz[Un] )
       !! ---------------------------------------------------------------------
-      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn 
+      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn
       INTEGER, DIMENSION(1:jpi          )     ::   ikbl
       REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
-      REAL(wp)                                ::   zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
-      REAL(wp)                                ::   zdU2,zdV2      
-      REAL(wp)                                ::   zwndi,zwndj
+      REAL(wp)                                ::   zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
+      REAL(wp)                                ::   zdU2, zdV2, zbuoy1, zbuoy2
+      REAL(wp)                                ::   zwndi, zwndj
       REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2
       REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF       
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF
       REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a
       REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b
@@ -648,5 +698,5 @@
       LOGICAL                                 ::   ln_Patankar    = .FALSE.
       LOGICAL                                 ::   ln_dumpvar     = .FALSE.
-      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl 
+      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl
       !
       tind  = nt_n
@@ -660,18 +710,18 @@
       !-------------
          !
-         ! Compute vertical shear         
+         ! Compute vertical shear
          DO jk = 2, jpkam1
-            DO ji = 1,jpi  
-               zcff        = 1.0_wp / e3w_abl( jk )**2                
-               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2                           
+            DO ji = 1,jpi
+               zcff        = 1.0_wp / e3w_abl( jk )**2
+               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
                zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
-               zsh2(ji,jk) = zdU+zdV
-            END DO
-         END DO   
+               zsh2(ji,jk) = zdU+zdV   !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 )
+            END DO
+         END DO
          !
          ! Compute brunt-vaisala frequency
          DO jk = 2, jpkam1
-            DO ji = 1,jpi 
-               zcff  = grav * itvref / e3w_abl( jk )  
+            DO ji = 1,jpi
+               zcff  = grav * itvref / e3w_abl( jk )
                zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta)
                zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        &
@@ -679,40 +729,45 @@
                zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi)
             END DO
-         END DO 
+         END DO
          !
          ! Terms for the tridiagonal problem
          DO jk = 2, jpkam1
-            DO ji = 1,jpi 
+            DO ji = 1,jpi
                zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point
-               zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation          
-               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
-               
+               zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
+               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
+
                z_elem_a( ji,     jk )   = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk   ) ! lower-diagonal
-               z_elem_c( ji,     jk )   = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal          
+               z_elem_c( ji,     jk )   = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal
                IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE
                   z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
-                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal       
-                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )             ! right-hand-side
+                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)      ! diagonal
+                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )               ! right-hand-side
                ELSE
                   z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
-                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal    
-                     &                     - e3w_abl(jk) * rDt_abl * zbuoy   
-                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )             ! right-hand-side                     
+                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)   &  ! diagonal
+                     &                     - e3w_abl(jk) * rDt_abl * zbuoy
+                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )             ! right-hand-side
                END IF
             END DO
-         END DO  
-                            
-         DO ji = 1,jpi    ! vector opt.            
-            zesrf   =  MAX( 4.63_wp * ustar2(ji,jj), tke_min )           
-            zetop   = tke_min             
-            z_elem_a ( ji,     1    ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp                                                   
-            z_elem_a ( ji,     jpka ) = 0._wp; z_elem_c ( ji,     jpka ) = 0._wp; z_elem_b ( ji,     jpka ) = 1._wp            
-            tke_abl( ji, jj,    1, nt_a ) = zesrf
-            tke_abl( ji, jj, jpka, nt_a ) = zetop 
-            zbn2(ji,jj,   1) = zbn2( ji,jj,   2) 
-            zsh2(ji,      1) = zsh2( ji,      2)                                      
-            zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 
+         END DO
+
+         DO ji = 1,jpi    ! vector opt.
+            zesrf   =  MAX( rn_esfc * ustar2(ji,jj), tke_min )
+            zetop   = tke_min
+
+            z_elem_a ( ji,        1 ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp
+            z_elem_a ( ji,     jpka ) = - 0.5 * rDt_abl * rn_Sch * (Avm_abl(ji,jj, jpka-1 )+Avm_abl(ji,jj, jpka ))  / e3t_abl( jpka )
+            z_elem_c ( ji,     jpka ) = 0._wp
+            z_elem_b ( ji,     jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka )
+
+            tke_abl  ( ji, jj,    1 , nt_a ) = zesrf
+            tke_abl  ( ji, jj, jpka , nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n )
+
+            zbn2(ji,jj,   1) = zbn2( ji,jj,   2)
+            zsh2(ji,      1) = zsh2( ji,      2)
+            zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1)
             zsh2(ji,   jpka) = zsh2( ji  , jpkam1)
-         END DO        
+         END DO
          !!
          !! Matrix inversion
@@ -720,9 +775,9 @@
          DO ji = 1,jpi
             zcff                  =  1._wp / z_elem_b( ji, 1 )
-            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       ) 
-            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a ) 
-         END DO             
-
-         DO jk = 2, jpka            
+            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       )
+            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a )
+         END DO
+
+         DO jk = 2, jpka
             DO ji = 1,jpi
                zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
@@ -732,12 +787,12 @@
             END DO
          END DO
-            
-         DO jk = jpkam1,1,-1            
+
+         DO jk = jpkam1,1,-1
             DO ji = 1,jpi
                tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a)
             END DO
          END DO
-         
-!!FL should not be needed because of Patankar procedure  
+
+!!FL should not be needed because of Patankar procedure
          tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
 
@@ -745,15 +800,15 @@
          !! Diagnose PBL height
          !! ----------------------------------------------------------
-         
-         
-         !                                                       
+
+
+         !
          ! arrays zRH, zFC and zCF are available at this point
          ! and zFC(:, 1 ) = 0.
          ! diagnose PBL height based on zsh2 and zbn2
          zFC (  :  ,1) = 0._wp
-         ikbl( 1:jpi ) = 0 
-         
+         ikbl( 1:jpi ) = 0
+
          DO jk = 2,jpka
-            DO ji = 1, jpi 
+            DO ji = 1, jpi
                zcff  = ghw_abl( jk-1 )
                zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
@@ -781,15 +836,15 @@
             ELSE
                pblh( ji, jj ) = ghw_abl(jpka)
-            END IF 
-         END DO
-      !-------------      
-      END DO      
-      !-------------
-      !
-      ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 
+            END IF
+         END DO
+      !-------------
+      END DO
+      !-------------
+      !
+      ! Optional : could add pblh smoothing if pblh is noisy horizontally ...
       IF(ln_smth_pblh) THEN
-         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
+         CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing)
          CALL smooth_pblh( pblh, msk_abl )
-         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)	  
+         CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing)
       ENDIF
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -799,92 +854,148 @@
       SELECT CASE ( nn_amxl )
       !
-      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom         
-#   define zlup zRH      
-#   define zldw zFC 
+      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom
+#   define zlup zRH
+#   define zldw zFC
          DO jj = 1, jpj     ! outer loop
             !
             DO ji = 1, jpi
-               mxl_abl ( ji, jj,    1 )  = 0._wp
-               mxl_abl ( ji, jj, jpka )  = mxl_min
-               zldw( ji,        1 )  = 0._wp              
-               zlup( ji,     jpka )  = 0._wp
-            END DO
-            !
-            DO jk = 2, jpkam1   
-               DO ji = 1, jpi  
-                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall )  
-                  mxl_abl( ji, jj, jk ) = MAX( mxl_min,  &
-                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   ) 
-               END DO
-            END DO   
+               mxld_abl ( ji, jj,    1 ) = 0._wp
+               mxld_abl ( ji, jj, jpka ) = mxl_min
+               mxlm_abl ( ji, jj,    1 ) = 0._wp
+               mxlm_abl ( ji, jj, jpka ) = mxl_min
+               zldw     ( ji,        1 ) = zrough(ji,jj) * rn_Lsfc
+               zlup     ( ji,     jpka ) = mxl_min
+            END DO
+            !
+            DO jk = 2, jpkam1
+               DO ji = 1, jpi
+                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall )
+                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min,  &
+                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   )
+               END DO
+            END DO
             !
             ! Limit mxl
-            DO jk = jpkam1,1,-1   
-               DO ji = 1, jpi 
-                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
-               END DO
-            END DO   
+            DO jk = jpkam1,1,-1
+               DO ji = 1, jpi
+                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
+               END DO
+            END DO
             !
             DO jk = 2, jpka
-               DO ji = 1, jpi  
-                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )   
-               END DO
-            END DO 
+               DO ji = 1, jpi
+                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
+               END DO
+            END DO
+            !
+!            DO jk = 1, jpka
+!               DO ji = 1, jpi
+!                  mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
+!                  mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ),  zlup( ji, jk ) )
+!               END DO
+!            END DO
             !
             DO jk = 1, jpka
                DO ji = 1, jpi
-                  mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )   
-               END DO
-            END DO                         
-            !
-         END DO
-#   undef zlup      
-#   undef zldw  
-         !
-         !
-      CASE ( 1 )             ! length-scale computed as the distance to the PBL height
-         DO jj = 1,jpj      ! outer loop
-            !
-            DO ji = 1, jpi   ! vector opt.
-               zcff = 1._wp / pblh( ji, jj )     ! inverse of hbl              
-               DO jk = 1, jpka              
-                  zsig  = MIN( zcff * ghw_abl( jk ), 1. )    
-                  zcff1 = pblh( ji, jj )                 
-                  mxl_abl( ji, jj, jk ) =  mxl_min                           &
-                     &  +  zsig         * (  amx1*zcff1 + bmx1*mxl_min ) &
-                     &  +  zsig *  zsig * (  amx2*zcff1 + bmx2*mxl_min ) &
-                     &  +  zsig**3      * (  amx3*zcff1 + bmx3*mxl_min ) &
-                     &  +  zsig**4      * (  amx4*zcff1 + bmx4*mxl_min ) &
-                     &  +  zsig**5      * (  amx5*zcff1 + bmx5*mxl_min )
-               END DO
-            END DO 
-            !            
-         END DO            
-         !
+                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
+!                  zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
+                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
+                  !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ),  zlup( ji, jk ) )
+                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
+               END DO
+            END DO
+            !
+         END DO
+#   undef zlup
+#   undef zldw
+         !
+         !
+      CASE ( 1 )           ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom
+#   define zlup zRH
+#   define zldw zFC
+         DO jj = 1, jpj     ! outer loop
+            !
+            DO jk = 2, jpkam1
+               DO ji = 1,jpi
+                              zcff        = 1.0_wp / e3w_abl( jk )**2
+                  zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
+                  zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
+                  zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
+                           ENDDO
+                        ENDDO
+                        !
+            DO ji = 1, jpi
+               zcff                      = zrough(ji,jj) * rn_Lsfc
+               mxld_abl ( ji, jj,    1 ) = zcff
+               mxld_abl ( ji, jj, jpka ) = mxl_min
+               mxlm_abl ( ji, jj,    1 ) = zcff
+               mxlm_abl ( ji, jj, jpka ) = mxl_min
+               zldw     ( ji,        1 ) = zcff
+               zlup     ( ji,     jpka ) = mxl_min
+            END DO
+            !
+            DO jk = 2, jpkam1
+               DO ji = 1, jpi
+                  zbuoy    = MAX( zbn2(ji, jj, jk), rsmall )
+                  zcff     = 2.*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) &
+                                &             + SQRT( rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.*zbuoy ) )
+                                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff )
+               END DO
+            END DO
+            !
+            ! Limit mxl
+            DO jk = jpkam1,1,-1
+               DO ji = 1, jpi
+                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
+               END DO
+            END DO
+            !
+            DO jk = 2, jpka
+               DO ji = 1, jpi
+                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
+               END DO
+            END DO
+            !
+            DO jk = 1, jpka
+               DO ji = 1, jpi
+                  !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
+                  !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
+                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
+                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
+                  !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )
+                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
+               END DO
+            END DO
+ !
+         END DO
+#   undef zlup
+#   undef zldw
+!         !
       CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale
          !
-#   define zlup zRH      
-#   define zldw zFC           
+#   define zlup zRH
+#   define zldw zFC
 ! zCF is used for matrix inversion
-!             
+!
        DO jj = 1, jpj      ! outer loop
-         
-         DO ji = 1, jpi           
-            zlup( ji,    1 ) = mxl_min
-            zldw( ji,    1 ) = mxl_min               
+
+         DO ji = 1, jpi
+            zcff             = zrough(ji,jj) * rn_Lsfc
+            zlup( ji,    1 ) = zcff
+            zldw( ji,    1 ) = zcff
             zlup( ji, jpka ) = mxl_min
-            zldw( ji, jpka ) = mxl_min                
-         END DO            
-         
+            zldw( ji, jpka ) = mxl_min
+         END DO
+
          DO jk = 2,jpka-1
             DO ji = 1, jpi
                zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
-               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)        
-            END DO
-         END DO         
+               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
+            END DO
+         END DO
          !!
          !! BL89 search for lup
-         !! ----------------------------------------------------------           
-         DO jk=2,jpka-1 
+         !! ----------------------------------------------------------
+         DO jk=2,jpka-1
             !
             DO ji = 1, jpi
@@ -892,17 +1003,20 @@
                zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
                ln_foundl(ji ) = .false.
-            END DO   
-            !           
+            END DO
+            !
             DO jkup=jk+1,jpka-1
                DO ji = 1, jpi
+                  zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall )
+                  zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall )
                   zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
-                     &               ( zbn2(ji,jj,jkup  )*(ghw_abl(jkup  )-ghw_abl(jk)) &
-                     &               + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
+                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk)) &
+                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )
                   IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
                      zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
-                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)                   
+                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
                      zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
-                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
-                     zlup(ji,jk) = zcff                
+                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )
+                     zlup(ji,jk) = zcff
+                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
                      ln_foundl(ji) = .true.
                   END IF
@@ -910,9 +1024,9 @@
             END DO
             !
-         END DO   
+         END DO
          !!
          !! BL89 search for ldwn
-         !! ----------------------------------------------------------          
-         DO jk=2,jpka-1         
+         !! ----------------------------------------------------------
+         DO jk=2,jpka-1
             !
             DO ji = 1, jpi
@@ -920,55 +1034,180 @@
                zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
                ln_foundl(ji ) = .false.
-            END DO  
-            !   
+            END DO
+            !
             DO jkdwn=jk-1,1,-1
-               DO ji = 1, jpi             
+               DO ji = 1, jpi
+                  zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
+                  zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
                   zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
-                     &               * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
-                                       + zbn2(ji,jj,jkdwn  )*(ghw_abl(jk)-ghw_abl(jkdwn  )) )   
-                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN 
+                     &               * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
+                                       + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )
+                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
                      zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
-                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )              
+                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )
                      zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
-                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
-                     zldw(ji,jk) = zcff 
-                     ln_foundl(ji) = .true.               
-                  END IF                                   
-               END DO
-            END DO
-            !      
+                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )
+                     zldw(ji,jk) = zcff
+                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
+                     ln_foundl(ji) = .true.
+                  END IF
+               END DO
+            END DO
+            !
          END DO
 
          DO jk = 1, jpka
-            DO ji = 1, jpi  
-               mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min )  
-            END DO
-         END DO  
+            DO ji = 1, jpi
+               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
+               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
+               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
+               mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
+            END DO
+         END DO
 
       END DO
-#   undef zlup      
-#   undef zldw           
-         !
-      END SELECT      
+#   undef zlup
+#   undef zldw
+         !
+     CASE ( 3 )           ! Bougeault & Lacarrere 89 length-scale
+         !
+#   define zlup zRH
+#   define zldw zFC
+! zCF is used for matrix inversion
+!
+       DO jj = 1, jpj      ! outer loop
+          !
+          DO jk = 2, jpkam1
+             DO ji = 1,jpi
+                            zcff        = 1.0_wp / e3w_abl( jk )**2
+                zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
+                zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
+                zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
+                         ENDDO
+                  ENDDO
+          zsh2(:,      1) = zsh2( :,      2)
+          zsh2(:,   jpka) = zsh2( :, jpkam1)
+
+                 DO ji = 1, jpi
+            zcff              = zrough(ji,jj) * rn_Lsfc
+                        zlup( ji,    1 )  = zcff
+            zldw( ji,    1 )  = zcff
+            zlup( ji, jpka ) = mxl_min
+            zldw( ji, jpka ) = mxl_min
+         END DO
+
+         DO jk = 2,jpka-1
+            DO ji = 1, jpi
+               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
+               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
+            END DO
+         END DO
+         !!
+         !! BL89 search for lup
+         !! ----------------------------------------------------------
+         DO jk=2,jpka-1
+            !
+            DO ji = 1, jpi
+               zCF(ji,1:jpka) = 0._wp
+               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
+               ln_foundl(ji ) = .false.
+            END DO
+            !
+            DO jkup=jk+1,jpka-1
+               DO ji = 1, jpi
+                  zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall )
+                  zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall )
+                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) *                 &
+                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk))      &
+                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )    &
+                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *        &
+                                         &               ( SQRT(tke_abl( ji, jj, jkup  , nt_a ))*zsh2(ji,jkup  ) &
+                                         &               + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) )
+
+                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
+                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
+                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
+                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
+                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )
+                     zlup(ji,jk) = zcff
+                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
+                     ln_foundl(ji) = .true.
+                  END IF
+               END DO
+            END DO
+            !
+         END DO
+         !!
+         !! BL89 search for ldwn
+         !! ----------------------------------------------------------
+         DO jk=2,jpka-1
+            !
+            DO ji = 1, jpi
+               zCF(ji,1:jpka) = 0._wp
+               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
+               ln_foundl(ji ) = .false.
+            END DO
+            !
+            DO jkdwn=jk-1,1,-1
+               DO ji = 1, jpi
+                  zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
+                  zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
+                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
+                     &               * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1))    &
+                     &                 +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )  &
+                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *          &
+                                         &               ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) &
+                                         &               + SQRT(tke_abl( ji, jj, jkdwn  , nt_a ))*zsh2(ji,jkdwn  ) )
+
+                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
+                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
+                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )
+                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
+                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )
+                     zldw(ji,jk) = zcff
+                                         zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
+                     ln_foundl(ji) = .true.
+                  END IF
+               END DO
+            END DO
+            !
+         END DO
+
+         DO jk = 1, jpka
+            DO ji = 1, jpi
+               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
+               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
+               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
+               mxld_abl( ji, jj, jk ) = MAX(  MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
+            END DO
+         END DO
+
+      END DO
+#   undef zlup
+#   undef zldw
+         !
+         !
+         !
+      END SELECT
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  Finalize the computation of turbulent visc./diff.
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      
+
       !-------------
       DO jj = 1, jpj     ! outer loop
       !-------------
-         DO jk = 1, jpka   
+         DO jk = 1, jpka
             DO ji = 1, jpi  ! vector opt.
-               zcff              = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk )  &
-                  &                                       * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) ) 
-               zcff2             =  1. / ( 1. + zcff )   !<-- phi_z(z)
-               zcff              = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
-!!FL: MAX function probably useless because of the definition of mxl_min             
+               zcff  = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk )  &
+               &     * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) )
+               zcff2 =  1. / ( 1. + zcff )   !<-- phi_z(z)
+               zcff  = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
+               !!FL: MAX function probably useless because of the definition of mxl_min
                Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
-               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )                              
-            END DO
-         END DO
-      !-------------         
-      END DO      
+               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )
+            END DO
+         END DO
+      !-------------
+      END DO
       !-------------
 
@@ -988,10 +1227,10 @@
       !!
       !! ---------------------------------------------------------------------
-	  REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk   
-	  REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
       INTEGER                                     :: ji,jj
-	  REAL(wp)                                    :: smth_a, smth_b
-	  REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY
-	  REAL(wp)                                    :: zumsk,zvmsk
+      REAL(wp)                                    :: smth_a, smth_b
+      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY
+      REAL(wp)                                    :: zumsk,zvmsk
       !!
       !!=========================================================
@@ -1005,15 +1244,15 @@
          zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk
       END_2D
-      
-	  DO_2D_10_11
+
+      DO_2D_10_11
          zvmsk = msk(ji,jj) * msk(ji,jj+1)
          zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk
-	  END_2D
-      
-	  DO_2D_10_00
+      END_2D
+
+      DO_2D_10_00
          zFY ( ji, jj  ) =   zdY ( ji, jj   )                        &
             & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   &
             &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  )
-	  END_2D
+      END_2D
 
       DO_2D_00_10
@@ -1029,5 +1268,5 @@
   &                 +zFY( ji, jj ) - zFY( ji, jj-1 )  )
       END_2D
-	  !!
+
 !---------------------------------------------------------------------------------------------------
    END SUBROUTINE smooth_pblh
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablrst.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablrst.F90	(revision 12587)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablrst.F90	(revision 12588)
@@ -117,5 +117,4 @@
       CALL iom_rstput( iter, nitrst, numraw, 'avm_abl', avm_abl(:,:,:           ) )
       CALL iom_rstput( iter, nitrst, numraw, 'avt_abl', avt_abl(:,:,:           ) )
-      CALL iom_rstput( iter, nitrst, numraw, 'mxl_abl', mxl_abl(:,:,:           ) )
       CALL iom_rstput( iter, nitrst, numraw,    'pblh',    pblh(:,:             ) )
       !
@@ -172,5 +171,4 @@
       CALL iom_get( numrar, jpdom_autoglo, 'avm_abl', avm_abl(:,:,:           ) )
       CALL iom_get( numrar, jpdom_autoglo, 'avt_abl', avt_abl(:,:,:           ) )
-      CALL iom_get( numrar, jpdom_autoglo, 'mxl_abl', mxl_abl(:,:,:           ) )
       CALL iom_get( numrar, jpdom_autoglo,    'pblh',    pblh(:,:             ) )
       CALL iom_delay_rst( 'READ', 'ABL', numrar )   ! read only abl delayed global communication variables
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/par_abl.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/par_abl.F90	(revision 12587)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/par_abl.F90	(revision 12588)
@@ -28,9 +28,11 @@
    LOGICAL , PUBLIC            ::   ln_hpgls_frc   !: forcing of ABL winds by large-scale pressure gradient 
    LOGICAL , PUBLIC            ::   ln_smth_pblh   !: smoothing of atmospheric PBL height 
+   LOGICAL , PUBLIC            ::   ln_topbc_neumann = .FALSE.  !: smoothing of atmospheric PBL height    
 
-   CHARACTER(len=256), PUBLIC ::   cn_ablrst_in     !: suffix of abl restart name (input)
-   CHARACTER(len=256), PUBLIC ::   cn_ablrst_out    !: suffix of abl restart name (output)
-   CHARACTER(len=256), PUBLIC ::   cn_ablrst_indir  !: abl restart input directory
-   CHARACTER(len=256), PUBLIC ::   cn_ablrst_outdir !: abl restart output directory
+   LOGICAL           , PUBLIC  ::   ln_rstart_abl    !: (de)activate abl restart
+   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_in     !: suffix of abl restart name (input)
+   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_out    !: suffix of abl restart name (output)
+   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_indir  !: abl restart input directory
+   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_outdir !: abl restart output directory
 
    !!---------------------------------------------------------------------
@@ -45,10 +47,12 @@
    REAL(wp), PUBLIC, PARAMETER ::   rn_Cek    = 258._wp                   !: Ekman constant for Richardson number 
    REAL(wp), PUBLIC, PARAMETER ::   rn_epssfc = 1._wp / ( 1._wp + 2.8_wp * 2.8_wp )
-   REAL(wp), PUBLIC            ::   rn_ceps                       !: namelist parameter
-   REAL(wp), PUBLIC            ::   rn_cm                         !: namelist parameter
-   REAL(wp), PUBLIC            ::   rn_ct                         !: namelist parameter
-   REAL(wp), PUBLIC            ::   rn_ce                         !: namelist parameter 
+   REAL(wp), PUBLIC            ::   rn_Ceps                       !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_Cm                         !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_Ct                         !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_Ce                         !: namelist parameter 
    REAL(wp), PUBLIC            ::   rn_Rod                        !: namelist parameter   
    REAL(wp), PUBLIC            ::   rn_Sch    
+   REAL(wp), PUBLIC            ::   rn_Esfc
+   REAL(wp), PUBLIC            ::   rn_Lsfc
    REAL(wp), PUBLIC            ::   mxl_min    
    REAL(wp), PUBLIC            ::   rn_ldyn_min                   !: namelist parameter
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90	(revision 12587)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90	(revision 12588)
@@ -68,15 +68,16 @@
       LOGICAL            ::   lluldl
       NAMELIST/namsbc_abl/ cn_dir, cn_dom, cn_ablrst_in, cn_ablrst_out,           &
-         &                 cn_ablrst_indir, cn_ablrst_outdir,                     &
+         &                 cn_ablrst_indir, cn_ablrst_outdir, ln_rstart_abl,      &
          &                 ln_hpgls_frc, ln_geos_winds, nn_dyn_restore,           &
          &                 rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max,   &
-         &                 nn_amxl, rn_cm, rn_ct, rn_ce, rn_ceps, rn_Rod, rn_Ric, &
+         &                 nn_amxl, rn_Cm, rn_Ct, rn_Ce, rn_Ceps, rn_Rod, rn_Ric, &
          &                 ln_smth_pblh
       !!---------------------------------------------------------------------
 
-      ! Namelist namsbc_abl in reference namelist : ABL parameters
+                                        ! Namelist namsbc_abl in reference namelist : ABL parameters
       READ  ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 )
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist' )
-      ! Namelist namsbc_abl in configuration namelist : ABL parameters
+      !
+                                        ! Namelist namsbc_abl in configuration namelist : ABL parameters
       READ  ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' )
@@ -165,5 +166,7 @@
       rn_Sch  = rn_ce / rn_cm
       mxl_min = (avm_bak / rn_cm) / sqrt( tke_min )
-
+      rn_Esfc =  1._wp / SQRT(rn_cm*rn_ceps)
+      rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm
+	  
       IF(lwp) THEN
          WRITE(numout,*)
@@ -171,5 +174,7 @@
          WRITE(numout,*) '    ~~~~~~~~~~~'
          IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale '
-         IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height '
+         IF(nn_amxl==1) WRITE(numout,*) 'Modified Deardorff 80 length-scale '
+         IF(nn_amxl==1) WRITE(numout,*) 'Bougeault and Lacarrere length-scale '		 
+         IF(nn_amxl==3) WRITE(numout,*) 'Rodier et al. length-scale '	
          WRITE(numout,*) ' Minimum value of atmospheric TKE           = ',tke_min,' m^2 s^-2'
          WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m'
@@ -178,4 +183,6 @@
          WRITE(numout,*) ' Constant for Schmidt number                = ',rn_Sch
          WRITE(numout,*) ' Constant for TKE dissipation               = ',rn_Ceps
+         WRITE(numout,*) ' Constant for TKE sfc boundary condition    = ',rn_Esfc
+         WRITE(numout,*) ' Constant for mxl sfc boundary condition    = ',rn_Lsfc
       END IF
 
@@ -248,7 +255,4 @@
          zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax
          rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8
-         !!GS: alternative shape
-         !rest_eq(:,:) = SIN( 0.5_wp*rpi*(zcff - ABS(ff_t(:,:))) / (zcff - 3.e-5) )**8
-         !WHERE(ABS(ff_t(:,:)).LE.3.e-5) rest_eq(:,:) = 1._wp
       ELSE
          rest_eq(:,:) = 1._wp
@@ -267,11 +271,11 @@
 
       ! initialize ABL from data or restart
-      IF( ln_rstart ) THEN
+      IF( ln_rstart_abl ) THEN
          CALL abl_rst_read
       ELSE
          CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step
 
-         u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:)
-         v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:)
+          u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:)
+          v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:)
          tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:)
          tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:)
@@ -280,5 +284,4 @@
          avm_abl(:,:,:          ) = avm_bak
          avt_abl(:,:,:          ) = avt_bak
-         mxl_abl(:,:,:          ) = mxl_min
          pblh   (:,:            ) = ghw_abl( 3 )  !<-- assume that the pbl contains 3 grid points
          u_abl  (:,:,:,nt_a     ) = 0._wp
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icesbc.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icesbc.F90	(revision 12587)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icesbc.F90	(revision 12588)
@@ -148,5 +148,5 @@
       CASE( jp_blk, jp_abl )  !--- bulk formulation & ABL formulation
                                   CALL blk_ice_2    ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),    &
-            &                                           sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )    ! 
+            &                                           sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )
          IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
Index: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90	(revision 12587)
+++ NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90	(revision 12588)
@@ -104,7 +104,7 @@
    REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements
    !
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme)
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme and ABL)
+   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice
+   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
 
    LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB
@@ -113,4 +113,5 @@
    LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
    LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB
+   LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature
    !
    INTEGER  ::   nhumi          ! choice of the bulk algorithm
@@ -170,5 +171,5 @@
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm
          &                 cn_dir , rn_zqt, rn_zu,                                    &
-         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
+         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, ln_tpot,  &
          &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB
       !!---------------------------------------------------------------------
@@ -593,7 +594,11 @@
          !#LB: because AGRIF hates functions that return something else than a scalar, need to
          !     use scalar version of gamma_moist() ...
-         DO_2D_11_11
-            ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt
-         END_2D
+         IF( ln_tpot ) THEN
+            DO_2D_11_11
+               ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt
+            END_2D
+         ELSE
+            ztpot = ptair(:,:)
+         ENDIF
       ENDIF
 
@@ -663,8 +668,8 @@
       ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
          CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
-            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         &
-            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 &
-            &               taum(:,:), psen(:,:), zqla(:,:),                  &
-            &               pEvap=pevp(:,:), prhoa=rhoa(:,:) )
+            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          &
+            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                  &
+            &               taum(:,:), psen(:,:), zqla(:,:),                   &
+            &               pEvap=pevp(:,:), prhoa=rhoa(:,:)                   )
 
          zqla(:,:) = zqla(:,:) * tmask(:,:,1)
@@ -893,5 +898,4 @@
 
       ! local scalars ( place there for vector optimisation purposes)
-      !IF (ln_abl) rhoa  (:,:)  = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI)
       zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
 
@@ -913,5 +917,5 @@
          IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   &
             &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' )
-      ELSE
+      ELSE ! ln_abl
          zztmp1 = 11637800.0_wp
          zztmp2 =    -5897.8_wp
