Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-ice.xml	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-ice.xml	(revision 11363)
@@ -90,4 +90,6 @@
           <field id="qemp_ice"     long_name="Downward Heat Flux from E-P over ice"                                                                              unit="W/m2" />
           <field id="albedo"       long_name="Mean albedo over sea ice and ocean"                                                                                unit=""     />
+          <field id="Cd_ice"       long_name="Momentum turbulent exchange coefficient"                                                                           unit=""     />
+          <field id="Ch_ice"       long_name="Heat turbulent exchange coefficient"                                                                               unit=""     />
 	  
 	  <!-- trends -->
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11363)
@@ -24,7 +24,7 @@
 
         <field id="toce"         long_name="temperature"                         standard_name="sea_water_potential_temperature"   unit="degC"     grid_ref="grid_T_3D"/>
-        <field id="toce_e3t"     long_name="temperature (thickness weighted)"                                                      unit="degC"     grid_ref="grid_T_3D" > toce * e3t </field >
+        <field id="toce_e3t"     long_name="temperature (thickness weighted)"                                                      unit="degC"     grid_ref="grid_T_3D" > toce * e3t </field>
         <field id="soce"         long_name="salinity"                            standard_name="sea_water_practical_salinity"      unit="1e-3"     grid_ref="grid_T_3D"/>
-        <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                      unit="1e-3"     grid_ref="grid_T_3D" > soce * e3t </field >
+        <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                      unit="1e-3"     grid_ref="grid_T_3D" > soce * e3t </field>
 
         <!-- t-eddy viscosity coefficients (ldfdyn) -->
@@ -33,5 +33,5 @@
 
         <field id="sst"          long_name="sea surface temperature"                            standard_name="sea_surface_temperature"             unit="degC"     />
-        <field id="sst2"         long_name="square of sea surface temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst * sst </field >
+        <field id="sst2"         long_name="square of sea surface temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst * sst </field>
         <field id="sstmax"       long_name="max of sea surface temperature"   field_ref="sst"   operation="maximum"                                                 />
         <field id="sstmin"       long_name="min of sea surface temperature"   field_ref="sst"   operation="minimum"                                                 />
@@ -45,5 +45,5 @@
         
         <field id="sss"          long_name="sea surface salinity"                               standard_name="sea_surface_salinity"                unit="1e-3"     />
-        <field id="sss2"         long_name="square of sea surface salinity"                                                                         unit="1e-6"      > sss * sss </field >
+        <field id="sss2"         long_name="square of sea surface salinity"                                                                         unit="1e-6"      > sss * sss </field>
         <field id="sssmax"       long_name="max of sea surface salinity"      field_ref="sss"   operation="maximum"                                                 />
         <field id="sssmin"       long_name="min of sea surface salinity"      field_ref="sss"   operation="minimum"                                                 />
@@ -54,5 +54,5 @@
 
         <field id="ssh"          long_name="sea surface height"                                 standard_name="sea_surface_height_above_geoid"             unit="m" />
-        <field id="ssh2"         long_name="square of sea surface height"                       standard_name="square_of_sea_surface_height_above_geoid"   unit="m2" > ssh * ssh </field >
+        <field id="ssh2"         long_name="square of sea surface height"                       standard_name="square_of_sea_surface_height_above_geoid"   unit="m2" > ssh * ssh </field>
         <field id="wetdep"       long_name="wet depth"                                          standard_name="wet_depth"                                  unit="m" />
         <field id="sshmax"       long_name="max of sea surface height"        field_ref="ssh"   operation="maximum"                                                 />
@@ -352,10 +352,13 @@
       <field_group id="ABL" > <!-- time step automaticaly defined based on nn_fsbc -->
 	<!-- variables available with ABL on atmospheric T grid-->
-	<field_group id="grid_ABL" grid_ref="grid_TA_3D" >
+	<field_group id="grid_ABL3D" grid_ref="grid_TA_3D" >
           <field id="u_abl"      long_name="ABL i-horizontal velocity"     standard_name="abl_x_velocity" unit="m/s"      />
           <field id="v_abl"      long_name="ABL j-horizontal velocity"     standard_name="abl_y_velocity" unit="m/s"      />
-          <field id="t_abl"      long_name="ABL potential temperature"  standard_name="abl_theta" unit="K"        />
-          <field id="q_abl"      long_name="ABL specific humidity"         standard_name="abl_qspe" unit="kg/kg"    />
-          <field id="pblh"       long_name="ABL height"                          standard_name="abl_height" unit="m"    grid_ref="grid_T_2D"     />
+          <field id="t_abl"      long_name="ABL potential temperature"     standard_name="abl_theta"      unit="K"        />
+          <field id="q_abl"      long_name="ABL specific humidity"         standard_name="abl_qspe"       unit="kg/kg"    />
+          <field id="coeft"      long_name="ABL nudging coefficient"       standard_name="coeft"          unit=""         />
+	</field_group>
+	<field_group id="grid_ABL2D" grid_ref="grid_TA_2D" >
+          <field id="pblh"       long_name="ABL height"                    standard_name="abl_height"     unit="m"        />
 	</field_group>
       </field_group> <!-- ABL -->
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/grid_def_nemo.xml
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/grid_def_nemo.xml	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/grid_def_nemo.xml	(revision 11363)
@@ -59,10 +59,10 @@
 
        <!-- ABL grid definition -->
+       <grid id="grid_TA_2D">
+         <domain domain_ref="grid_T" />
+       </grid>
        <grid id="grid_TA_3D">
          <domain domain_ref="grid_T" />
          <axis  id="ght_abl" />
-       </grid>
-       <grid id="grid_TA_2D">
-         <domain domain_ref="grid_T" />
        </grid>
        <grid id="grid_WA_3D">
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref	(revision 11363)
@@ -283,5 +283,5 @@
 &namsbc_abl    !   Atmospheric Boundary Layer formulation           (ln_abl = T)
 !-----------------------------------------------------------------------
-   cn_dir         = 'atm_forcing/'      !  root directory for the location of the ABL grid file
+   cn_dir         = './'      !  root directory for the location of the ABL grid file
    cn_dom         = 'dom_cfg_abl.nc'
    ln_hpgls_frc   = .false.
@@ -301,5 +301,5 @@
    rn_Rod        = 0.15       ! c0 in RMCA17 mixing length formulation (not yet implemented)
    rn_Ric        = 0.139      !  Critical Richardson number (to compute PBL height and diffusivities)
-   ln_smth_pblh  = .true.    !  Smoothing of PBL height with a 2d Hanning filter
+   ln_smth_pblh  = .false.   !  Smoothing of PBL height with a 2d Hanning filter
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11363)
@@ -153,6 +153,6 @@
          DO jk = 3, jpkam1
             DO ji = 1, jpi   ! vector opt.
-               z_elem_a( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )          ! lower-diagonal
-               z_elem_c( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )          ! upper-diagonal       
+               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
@@ -162,7 +162,7 @@
             ! Neumann at the bottom          
             z_elem_a( ji,     2              ) = 0._wp
-            z_elem_c( ji,     2              ) = - rdt * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                         
+            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 * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
+            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 )
@@ -185,6 +185,6 @@
                   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 * zztmp1				  
-                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2               
+				  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 
@@ -197,6 +197,6 @@
 				  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 * zztmp1
-                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2               
+                  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                                     
@@ -244,14 +244,14 @@
          DO jj = 1, jpj
             DO ji = 1, jpi           
-               zcff = ( fft_abl(ji,jj) * rdt )*( fft_abl(ji,jj) * rdt )  ! (f dt)**2
+               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 * fft_abl(ji, jj) * v_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 * fft_abl(ji, jj) * u_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)                
             END DO
@@ -267,7 +267,7 @@
                DO ji = 1, jpi 
                   u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
-                     &                      - rdt * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
+                     &                      - 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 * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
+                     &                      + rdt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
                END DO
             END DO
@@ -280,6 +280,6 @@
             DO jk = 1, jpka
                DO ji = 1, jpi 
-                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgu_dta(ji,jj,jk)
-                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgv_dta(ji,jj,jk)  
+                  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)  
                ENDDO
             ENDDO
@@ -298,6 +298,6 @@
          DO jk = 3, jpkam1
             DO ji = 1, jpi  
-               z_elem_a( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal
-               z_elem_c( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal                
+               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_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal
             END DO
@@ -307,5 +307,5 @@
             !++ Surface boundary condition
             z_elem_a( ji,     2    ) = 0._wp
-            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
+            z_elem_c( ji,     2    ) = - rdt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
             !
 			zztmp1  = pcd_du(ji, jj)
@@ -316,9 +316,9 @@
 			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 * zztmp1 			
-            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rdt * zztmp2
+			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 * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka ) 
+            !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 )                                               
@@ -365,6 +365,6 @@
          DO jk = 3, jpkam1
             DO ji = 1, jpi   
-               z_elem_a( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
-               z_elem_c( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal              
+               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_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal
             END DO
@@ -374,5 +374,5 @@
             !++ Surface boundary condition
             z_elem_a( ji,     2    ) = 0._wp
-            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )        
+            z_elem_c( ji,     2    ) = - rdt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )        
             !
 			zztmp1 = pcd_du(ji, jj)
@@ -383,8 +383,8 @@
 			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 * zztmp1 
-            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rdt * zztmp2
+            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 * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
+            !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 )                                               
@@ -434,13 +434,13 @@
             DO jj = 2, jpj
                DO ji = 2, jpi
-                  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) ) 
+                  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) ) 
                   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 * rdt * zcff2   ! zcff = 1 for masked points
-                  
-				  zcff  = zcff * rest_eq(ji,jj) ; z_cft( ji, jj, jk ) = zcff
+                  zcff  = (1._wp-zmsk) + zmsk * zcff2 * rdt   ! zcff = 1 for masked points
+                                                              ! rdt = 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 )   &
@@ -448,5 +448,5 @@
                   v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   &
                      &                               + zcff   * pv_dta( ji, jj, jk       )
-                END DO
+               END DO
             END DO   
          !-------------
@@ -460,11 +460,12 @@
          DO jj = 1,jpj
             DO ji = 1,jpi 
-               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) ) 
-               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 * rdt * zcff2   ! zcff = 1 for masked points      
+               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) ) 
+               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 * rdt   ! zcff = 1 for masked points
+                                                           ! rdt = 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 )
@@ -483,6 +484,6 @@
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 	  !
-      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...
+      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...
 	  !
       CALL iom_put ( 'u_abl',  u_abl(:,:,2:jpka,nt_a      ))
@@ -490,5 +491,5 @@
       CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta))
       CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa))
-!      CALL iom_put ( 'coeft',  z_cft(:,:,2:jpka           ))
+      CALL iom_put ( 'coeft',  z_cft(:,:,2:jpka           ))
       CALL iom_put (  'pblh',   pblh(:,:                  ))
       !
@@ -676,15 +677,15 @@
                zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
                
-               z_elem_a( ji,     jk )   = - 0.5_wp * rdt * 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 * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal          
+               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          
                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 * 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 * ( zbuoy + zshear ) )             ! right-hand-side
+                     &                     + 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
                ELSE
                   z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
-                     &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal    
-                     &                     - e3w_abl(jk) * rdt * zbuoy   
-                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt *  zshear )             ! right-hand-side                     
+                     &                     + 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                     
                END IF
             END DO
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90	(revision 11363)
@@ -19,4 +19,5 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_qa  = 2     !: indice for humidity
    INTEGER , PUBLIC, PARAMETER ::   jptime = 2     !: number of time indices stored in memory
+
    !!---------------------------------------------------------------------
    !! namABL Namelist options  
@@ -69,4 +70,6 @@
    ! parameter for the semi-implicit treatment of Coriolis term  
    REAL(wp), PUBLIC, PARAMETER ::   gamma_Cor  = 0.55_wp
+   ! ABL timestep
+   REAL(wp), PUBLIC            :: rdt_abl
 
    !!---------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11363)
@@ -200,19 +200,24 @@
       jp_pblh_max = ghw_abl(jpka-3) / jp_bmax  !<-- at least 3 grid points at the top    have value rn_ltra_max
 
+      ! ABL timestep
+      rdt_abl = nn_fsbc * rdt
+
       ! Check parameters for dynamics
-      zcff  = jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2   &
-         &  + jp_alp1_dyn * jp_bmin    + jp_alp0_dyn 
-      zcff1 = jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2   &
-         &  + jp_alp1_dyn * jp_bmax    + jp_alp0_dyn      
+      zcff  = ( jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2   &
+         &    + jp_alp1_dyn * jp_bmin    + jp_alp0_dyn ) * rdt_abl
+      zcff1 = ( jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2   &
+         &    + jp_alp1_dyn * jp_bmax    + jp_alp0_dyn ) * rdt_abl
       IF(lwp) THEN
          IF(nn_dyn_restore > 0) THEN
-            WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff  * rdt
-            WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 * rdt
+            WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff
+            WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1
             ! Check that restoring coefficients are between 0 and 1
-            IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   &
+            !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   &
+            IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   &
                &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' )
-            IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   &
+            !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   &
+            IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   &
                &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' )
-            IF( zcff  * rdt > zcff1 * rdt ) &
+            IF( zcff  > zcff1                    )   &
                &                   CALL ctl_stop( 'abl_init : rn_ldyn_max should be smaller than rn_ldyn_min' )
          END IF
@@ -220,17 +225,19 @@
 
       ! Check parameters for active tracers
-      zcff  = jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2   &
-         &  + jp_alp1_tra * jp_bmin    + jp_alp0_tra 
-      zcff1 = jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   &
-         &  + jp_alp1_tra * jp_bmax    + jp_alp0_tra  
+      zcff  = ( jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2   &
+         &    + jp_alp1_tra * jp_bmin    + jp_alp0_tra ) * rdt_abl
+      zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   &
+         &    + jp_alp1_tra * jp_bmax    + jp_alp0_tra ) * rdt_abl 
       IF(lwp) THEN
-         WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff  * rdt
-         WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 * rdt
+         WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff
+         WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1
          ! Check that restoring coefficients are between 0 and 1
-         IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   &
+         !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   &
+         IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   &
             &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' )
-         IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   &
+         !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   &
+         IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   &
             &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' )
-         IF( zcff  * rdt > zcff1  * rdt ) &
+         IF( zcff  > zcff1                    )   &
             &                   CALL ctl_stop( 'abl_init : rn_ltra_max should be smaller than rn_ltra_min' )
       END IF
@@ -244,5 +251,8 @@
       IF( nn_dyn_restore == 1 ) THEN
          zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax 
-         rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff)/zcff ) )**8
+         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
Index: /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90	(revision 11362)
+++ /NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90	(revision 11363)
@@ -654,8 +654,7 @@
          CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!
+         !
          IF( ln_abl ) THEN
-         !                                                                                      !!! nid_A : 3D
-			CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl
+            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl
                &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
             CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl
@@ -677,8 +676,7 @@
                &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )
 #endif
-		    CALL histend( nid_A, snc4chunks=snc4set )
-		 !
-		 ENDIF
-!
+            CALL histend( nid_A, snc4chunks=snc4set )
+         ENDIF
+         !
          IF( ln_icebergs ) THEN
             CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
@@ -838,28 +836,28 @@
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
-!
+      !
       IF( ln_abl ) THEN 
-	     ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
-	     IF( ln_mskland )   THEN 
-		    DO jk=1,jpka
-		       zw3d_abl(:,:,jk) = tmask(:,:,1)
-            END DO
-		 ELSE
+         ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
+         IF( ln_mskland )   THEN 
+            DO jk=1,jpka
+               zw3d_abl(:,:,jk) = tmask(:,:,1)
+            END DO       
+         ELSE
             zw3d_abl(:,:,:) = 1._wp		 
          ENDIF			
-		 CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh 
-	     CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
-	     CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
-	     CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
-	     CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl		 
-	     CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
-	     CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
-	     CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl	
+         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh 
+         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
+         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
+         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
+         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl		 
+         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
+         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
+         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl	
 #if defined key_si3
          CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i
 #endif
-		 DEALLOCATE(zw3d_abl)
-	  ENDIF
-!
+         DEALLOCATE(zw3d_abl)
+      ENDIF
+      !
       IF( ln_icebergs ) THEN
          !
