Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90	(revision 11273)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90	(revision 11275)
@@ -29,7 +29,7 @@
    USE lib_mpp           ! MPP library
 #if defined key_iomput
-   USE sbc_oce  , ONLY :   nn_fsbc         ! ocean space and time domain
-   USE trc_oce  , ONLY :   nn_dttrc        !  !: frequency of step on passive tracers
-   USE icb_oce  , ONLY :   nclasses, class_num       !  !: iceberg classes
+   USE sbc_oce  , ONLY :   nn_fsbc, ght_abl, ghw_abl, e3w_abl, jpka, jpkam1
+   USE trc_oce  , ONLY :   nn_dttrc                  ! frequency of step on passive tracers
+   USE icb_oce  , ONLY :   nclasses, class_num       ! iceberg classes
 #if defined key_si3
    USE ice      , ONLY :   jpl
@@ -112,4 +112,5 @@
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zt_bnds, zw_bnds
+      REAL(wp), DIMENSION(2   ,jpkam1)      ::   za_bnds     ! ABL vertical boundaries
       LOGICAL ::   ll_tmppatch = .TRUE.    !: seb: patch before we remove periodicity
       INTEGER ::   nldi_save, nlei_save    !:      and close boundaries in output files
@@ -195,9 +196,18 @@
       ! vertical grid definition
       IF(.NOT.llrst_context) THEN
-          CALL iom_set_axis_attr( "deptht",  paxis = gdept_1d )
-          CALL iom_set_axis_attr( "depthu",  paxis = gdept_1d )
-          CALL iom_set_axis_attr( "depthv",  paxis = gdept_1d )
-          CALL iom_set_axis_attr( "depthw",  paxis = gdepw_1d )
-
+          CALL iom_set_axis_attr(  "deptht", paxis = gdept_1d )
+          CALL iom_set_axis_attr(  "depthu", paxis = gdept_1d )
+          CALL iom_set_axis_attr(  "depthv", paxis = gdept_1d )
+          CALL iom_set_axis_attr(  "depthw", paxis = gdepw_1d )
+
+          ! ABL
+          IF( .NOT. ALLOCATED(ght_abl) ) THEN   ! force definition for xml files (xios) 
+             ALLOCATE( ght_abl(jpka), ghw_abl(jpka), e3t_abl(jpka), e3w_abl(jpka) )   ! default allocation needed by iom
+             ght_abl(:) = -1._wp   ;   ghw_abl(:) = -1._wp
+             e3t_abl(:) = -1._wp   ;   e3w_abl(:) = -1._wp
+          ENDIF
+          CALL iom_set_axis_attr( "ght_abl", ght_abl(2:jpka) )
+          CALL iom_set_axis_attr( "ghw_abl", ghw_abl(2:jpka) )
+          
           ! Add vertical grid bounds
           jkmin = MIN(2,jpk)  ! in case jpk=1 (i.e. sas2D)
@@ -208,8 +218,16 @@
           zw_bnds(2,1:jpkm1  ) = gdepw_1d(jkmin:jpk)
           zw_bnds(2,jpk:     ) = gdepw_1d(jpk) + e3t_1d(jpk)
-          CALL iom_set_axis_attr( "deptht", bounds=zw_bnds )
-          CALL iom_set_axis_attr( "depthu", bounds=zw_bnds )
-          CALL iom_set_axis_attr( "depthv", bounds=zw_bnds )
-          CALL iom_set_axis_attr( "depthw", bounds=zt_bnds )
+          CALL iom_set_axis_attr(  "deptht", bounds=zw_bnds )
+          CALL iom_set_axis_attr(  "depthu", bounds=zw_bnds )
+          CALL iom_set_axis_attr(  "depthv", bounds=zw_bnds )
+          CALL iom_set_axis_attr(  "depthw", bounds=zt_bnds )
+
+          ! ABL
+          za_bnds(1,:) = ghw_abl(1:jpkam1)
+          za_bnds(2,:) = ghw_abl(2:jpka  )
+          CALL iom_set_axis_attr( "ght_abl", bounds=za_bnds )
+          za_bnds(1,:) = ght_abl(2:jpka  )
+          za_bnds(2,:) = ght_abl(2:jpka  ) + e3w_abl(2:jpka)
+          CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds )
           !
 # if defined key_floats
@@ -1943,12 +1961,10 @@
       !
       INTEGER :: ji, jj, jn, ni, nj
-      INTEGER :: icnr, jcnr                                    ! Offset such that the vertex coordinate (i+icnr,j+jcnr)
-      !                                                        ! represents the bottom-left corner of cell (i,j)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j)
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_fld       ! Working array to determine where to rotate cells
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_rot       ! Lat/lon working array for rotation of cells
-      !!----------------------------------------------------------------------
-      !
-      ALLOCATE( z_bnds(4,jpi,jpj,2), z_fld(jpi,jpj), z_rot(4,2)  )
+      INTEGER :: icnr, jcnr                             ! Offset such that the vertex coordinate (i+icnr,j+jcnr)
+      !                                                 ! represents the bottom-left corner of cell (i,j)
+      REAL(wp), DIMENSION(4,jpi,jpj,2) ::   z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j)
+      REAL(wp), DIMENSION(  jpi,jpj  ) ::   z_fld       ! Working array to determine where to rotate cells
+      REAL(wp), DIMENSION(4        ,2) ::   z_rot       ! Lat/lon working array for rotation of cells
+      !!----------------------------------------------------------------------
       !
       ! Offset of coordinate representing bottom-left corner
@@ -2031,6 +2047,4 @@
           &                                    bounds_lon = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,2),(/ 4,ni*nj /)), nvertex=4 )
       !
-      DEALLOCATE( z_bnds, z_fld, z_rot ) 
-      !
    END SUBROUTINE set_grid_bounds
 
@@ -2115,4 +2129,5 @@
       f_op%timestep = nn_fsbc  ;  f_of%timestep =  0  ; CALL iom_set_field_attr('SBC'             , freq_op=f_op, freq_offset=f_of)
       f_op%timestep = nn_fsbc  ;  f_of%timestep =  0  ; CALL iom_set_field_attr('SBC_scalar'      , freq_op=f_op, freq_offset=f_of)
+      f_op%timestep = nn_fsbc  ;  f_of%timestep =  0  ; CALL iom_set_field_attr('ABL'             , freq_op=f_op, freq_offset=f_of)
       f_op%timestep = nn_dttrc ;  f_of%timestep =  0  ; CALL iom_set_field_attr('ptrc_T'          , freq_op=f_op, freq_offset=f_of)
       f_op%timestep = nn_dttrc ;  f_of%timestep =  0  ; CALL iom_set_field_attr('diad_T'          , freq_op=f_op, freq_offset=f_of)
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbc_oce.F90	(revision 11273)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbc_oce.F90	(revision 11275)
@@ -11,4 +11,5 @@
    !!            4.0  ! 2012-05  (C. Rousset) add attenuation coef for use in ice model 
    !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk)
+   !!            4.0  ! 2019-03  (F. Lemarié, G. Samson) add compatibility with ABL mode    
    !!----------------------------------------------------------------------
 
@@ -34,4 +35,5 @@
    LOGICAL , PUBLIC ::   ln_flx         !: flux      formulation
    LOGICAL , PUBLIC ::   ln_blk         !: bulk formulation
+   LOGICAL , PUBLIC ::   ln_abl         !: Atmospheric boundary layer model
 #if defined key_oasis3
    LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used
@@ -77,6 +79,7 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_flx     = 2        !: flux                          formulation
    INTEGER , PUBLIC, PARAMETER ::   jp_blk     = 3        !: bulk                          formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 4        !: Pure ocean-atmosphere Coupled formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_none    = 5        !: for OPA when doing coupling via SAS module
+   INTEGER , PUBLIC, PARAMETER ::   jp_abl     = 4        !: Atmospheric boundary layer    formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 5        !: Pure ocean-atmosphere Coupled formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_none    = 6        !: for OPA when doing coupling via SAS module
    
    !!----------------------------------------------------------------------
@@ -107,10 +110,9 @@
    INTEGER , PUBLIC ::  ncpl_qsr_freq            !: qsr coupling frequency per days from atmosphere
    !
-   LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress)
    !!                                   !!   now    ! before   !!
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
-   !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads
+   !! wndm is used compute surface gases exchanges in ice-free ocean or leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr               !: sea heat flux:     solar                     [W/m2]
@@ -136,4 +138,12 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
+
+   !!---------------------------------------------------------------------
+   !! ABL Vertical Domain size  
+   !!---------------------------------------------------------------------
+   INTEGER , PUBLIC            ::   jpka   = 2     !: ABL number of vertical levels (default definition)
+   INTEGER , PUBLIC            ::   jpkam1 = 1     !: jpka-1
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   ght_abl, ghw_abl          !: ABL geopotential height (needed for iom)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_abl, e3w_abl          !: ABL vertical scale factors (needed for iom)
 
    !!----------------------------------------------------------------------
@@ -182,7 +192,7 @@
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      &
          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
-         !
+      !
       ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
-         !
+      !
       sbc_oce_alloc = MAXVAL( ierr )
       CALL mpp_sum ( 'sbc_oce', sbc_oce_alloc )
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcabl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcabl.F90	(revision 11275)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcabl.F90	(revision 11275)
@@ -0,0 +1,50 @@
+MODULE sbcabl
+   !!======================================================================
+   !!                       ***  MODULE  sbcabl  ***
+   !! Ocean forcing:  momentum, heat and freshwater flux formulation 
+   !!                 derived from an ABL model 
+   !!=====================================================================
+   !! History :  4.0  !  2019-03  (F. Lemarié & G. Samson)  Original code
+   !!----------------------------------------------------------------------
+   USE sbc_oce, ONLY :   ght_abl, ghw_abl, e3t_abl, e3w_abl
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_abl_init       ! routine called in sbcmod module
+   PUBLIC   sbc_abl            ! routine called in sbcmod module
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO-consortium (2014)
+   !! $Id: sbcabl.F90 6416 2016-04-01 12:22:17Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_abl_init
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_abl_init  ***
+      !!
+      !! ** Purposes :   dummy routine for compilation 
+      !!
+      !!----------------------------------------------------------------------
+      
+   END SUBROUTINE sbc_abl_init
+
+   
+   SUBROUTINE sbc_abl( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_abl  ***
+      !!
+      !! ** Purposes :   dummy routine for compilation 
+      !!
+      !!---------------------------------------------------------------------
+      INTEGER ,         INTENT(in) ::   kt   ! ocean time step
+      !!
+   END SUBROUTINE sbc_abl
+
+
+   !!======================================================================
+END MODULE sbcabl
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11273)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11275)
@@ -18,4 +18,5 @@
    !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
    !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle) 
+   !!            4.0  !  2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE)
    !!----------------------------------------------------------------------
 
@@ -23,5 +24,6 @@
    !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
    !!   sbc_blk       : bulk formulation as ocean surface boundary condition
-   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
+   !!       blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model  (ln_abl=TRUE)  
+   !!       blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step  (ln_abl=TRUE) 
    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
@@ -65,4 +67,8 @@
    PUBLIC   sbc_blk_init  ! called in sbcmod
    PUBLIC   sbc_blk       ! called in sbcmod
+   PUBLIC   blk_oce_1     ! called in sbcabl
+   PUBLIC   blk_oce_2     ! called in sbcabl
+   PUBLIC   rho_air       ! called in ablmod
+   PUBLIC   cp_air        ! called in ablmod
 #if defined key_si3
    PUBLIC   blk_ice_tau   ! routine called in icesbc
@@ -71,5 +77,9 @@
 #endif 
 
-!!Lolo: should ultimately be moved in the module with all physical constants ?
+  INTERFACE cp_air
+    MODULE PROCEDURE cp_air_0d, cp_air_2d
+  END INTERFACE
+
+  !!Lolo: should ultimately be moved in the module with all physical constants ?
 !!gm  : In principle, yes.
    REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg]
@@ -80,17 +90,18 @@
    REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
 
-   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
-   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin)
-   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
-   INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2)
-   INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2)
-   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
-   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
-   INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa)
-   INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point
-
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
+   INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read
+   INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   !: index of 10m wind velocity (i-component) (m/s)    at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   !: index of 10m wind velocity (j-component) (m/s)    at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   !: index of 10m air temperature             (Kelvin)
+   INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   !: index of specific humidity               ( % )
+   INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   !: index of solar heat                      (W/m2)
+   INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   !: index of Long wave                       (W/m2)
+   INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   !: index of total precipitation (rain+snow) (Kg/m2/s)
+   INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   !: index of snow (solid prcipitation)       (kg/m2/s)
+   INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   !: index of sea level pressure              (Pa)
+   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   !: index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
+   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
+
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   !  structure of input fields (file informations, fields read)
 
    !                                             !!! Bulk parameters
@@ -107,5 +118,4 @@
    LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
    !
-   LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data
    REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation
    REAL(wp) ::   rn_efac        ! multiplication factor for evaporation
@@ -161,5 +171,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
+      INTEGER  ::   ifpr                  ! dummy loop indice and argument
       INTEGER  ::   ios, ierror, ioptio   ! Local integer
       !!
@@ -168,9 +178,9 @@
       TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
-      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        "
+      TYPE(FLD_N) ::   sn_slp , sn_hpgi, sn_hpgj               !       "                        "
       NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
-         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
+         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       &
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
-         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
+         &                 cn_dir , rn_zqt, rn_zu,                                    & 
          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
       !!---------------------------------------------------------------------
@@ -214,22 +224,30 @@
       slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
       slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
-      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
-      !
-      lhftau = ln_taudif                     !- add an extra field if HF stress is used
-      jfld = jpfld - COUNT( (/.NOT.lhftau/) )
+      slf_i(jp_slp)  = sn_slp
+      slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj  
       !
       !                                      !- allocate the bulk structure
-      ALLOCATE( sf(jfld), STAT=ierror )
+      ALLOCATE( sf(jpfld), STAT=ierror )
       IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
-      DO ifpr= 1, jfld
-         ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
-         IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
-         IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   &
-            &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
-            &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
-
-      END DO
       !                                      !- fill the bulk structure with namelist informations
       CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
+      !
+      DO ifpr= 1, jpfld
+         !
+         IF( TRIM(sf(ifpr)%clrootname) /= 'NOT USED' ) THEN
+            IF( ln_abl .AND. (     ifpr == jp_wndi .OR. ifpr == jp_wndj .OR. ifpr == jp_humi   &
+               &              .OR. ifpr == jp_hpgi .OR. ifpr == jp_hpgj .OR. ifpr == jp_tair ) ) THEN
+               ALLOCATE( sf(ifpr)%fnow(jpi,jpj,jpka) )
+               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,jpka,2) )
+            ELSE
+               ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
+               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
+            ENDIF
+            !
+            IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   &
+               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
+               &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
+         ENDIF
+      END DO
       !
       IF ( ln_wave ) THEN
@@ -252,4 +270,10 @@
       ENDIF 
       !
+      IF( ln_abl ) THEN      ! ABL: read 3D fields for wind, temperature and humidity     
+         rn_zqt = ght_abl(2)          ! set the bulk altitude to ABL first level
+         rn_zu  = ght_abl(2)
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
+      ENDIF
       !           
       IF(lwp) THEN                     !** Control print
@@ -261,5 +285,4 @@
          WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0
          WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
-         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
          WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
          WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
@@ -294,12 +317,12 @@
       !!      the 10m wind velocity (i-component) (m/s)    at T-point
       !!      the 10m wind velocity (j-component) (m/s)    at T-point
-      !!      the 10m or 2m specific humidity     ( % )
+      !!      the 10m or 2m specific humidity     (kg/kg)
       !!      the solar heat                      (W/m2)
       !!      the Long wave                       (W/m2)
       !!      the 10m or 2m air temperature       (Kelvin)
       !!      the total precipitation (rain+snow) (Kg/m2/s)
-      !!      the snow (solid prcipitation)       (kg/m2/s)
-      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
-      !!              (2) CALL blk_oce
+      !!      the snow (solid precipitation)      (kg/m2/s)
+      !!      ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 
+      !!              (2) CALL blk_oce_1 and blk_oce_2
       !!
       !!      C A U T I O N : never mask the surface stress fields
@@ -318,4 +341,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp
       !!---------------------------------------------------------------------
       !
@@ -323,6 +348,18 @@
       !
       !                                            ! compute the surface ocean fluxes using bulk formulea
-      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
-
+      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+         !
+         CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in
+            &                zssq, zcd_du, zsen, zevp )                              !   =>> out
+         
+         CALL blk_oce_2( kt, sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in
+            &                sf(jp_snow)%fnow(:,:,1), sst_m,                     &   !   <<= in
+            &                zsen, zevp )                                            !   <=> in out
+         
+      ENDIF
+         
 #if defined key_cice
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
@@ -343,13 +380,179 @@
 
 
-   SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE blk_oce  ***
-      !!
-      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
-      !!      the ocean surface at each time step
-      !!
-      !! ** Method  :   bulk formulea for the ocean using atmospheric
-      !!      fields read in sbc_read
+   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, &  ! inp
+              &              pslp,  pst  , pu   , pv,    &  ! inp
+              &              pssq, pcd_du, psen, pevp    )  ! out
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_oce_1  ***
+      !!
+      !! ** Purpose :   Computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration
+      !!
+      !! ** Method  :   bulk formulae using atmospheric
+      !!                fields from the ABL model at previous time-step
+      !!
+      !! ** Outputs : - pssq    : surface humidity used to compute latent heat flux (kg/kg)
+      !!              - pcd_du  : Cd x |dU| at T-points  (m/s)
+      !!              - psen    : Ch x |dU| at T-points  (m/s)
+      !!              - pevp    : Ce x |dU| at T-points  (m/s)
+      !!---------------------------------------------------------------------
+      INTEGER , INTENT(in   )                 ::   kt     ! time step index
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndi  ! atmospheric wind at U-point              [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndj  ! atmospheric wind at V-point              [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   phumi  ! specific humidity at T-points            [kg/kg]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   ptair  ! potential temperature at T-points        [Kelvin]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pslp   ! sea-level pressure                       [Pa]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celcius]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg]
+      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   psen   ! Ch x |dU| at T-points                    [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pevp   ! Ce x |dU| at T-points                    [m/s]
+      !
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      REAL(wp) ::   zztmp                ! local variable
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
+      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin
+      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
+      !!---------------------------------------------------------------------
+      !
+      ! local scalars ( place there for vector optimisation purposes)
+      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
+
+      ! ----------------------------------------------------------------------------- !
+      !      0   Wind components and module at T-point relative to the moving ocean   !
+      ! ----------------------------------------------------------------------------- !
+
+      ! ... components ( U10m - U_oce ) at T-point (unmasked)
+#if defined key_cyclone
+      zwnd_i(:,:) = 0._wp
+      zwnd_j(:,:) = 0._wp
+      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vect. opt.
+            pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
+            pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
+         END DO
+      END DO
+#endif
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vect. opt.
+            zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
+            zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
+      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
+      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
+         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
+
+      ! ----------------------------------------------------------------------------- !
+      !     I    Turbulent FLUXES                                                    !
+      ! ----------------------------------------------------------------------------- !
+
+      ! ... specific humidity at SST and IST tmask(
+      pssq(:,:) = 0.98 * q_sat( zst(:,:), pslp(:,:) )
+      !
+      IF( ln_abl ) THEN
+         ztpot = ptair(:,:)
+      ELSE
+         ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
+         !    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
+         !    (since reanalysis products provide T at z, not theta !)
+         ztpot = ptair(:,:) + gamma_moist( ptair(:,:), phumi(:,:) ) * rn_zqt
+      ENDIF
+
+      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
+      !
+      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! NCAR-COREv2
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.0
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.5
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! ECMWF
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+      CASE DEFAULT
+         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
+      END SELECT
+
+
+      IF ( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
+!! FL do we need this multiplication by tmask ... ???
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zztmp = zU_zu(ji,jj) * tmask(ji,jj,1)
+               wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod 
+               pcd_du(ji,jj) = zztmp * Cd_atm(ji,jj)
+               psen(ji,jj)   = zztmp * Ch_atm(ji,jj)
+               pevp(ji,jj)   = zztmp * Ce_atm(ji,jj)       
+            END DO
+         END DO
+
+      ELSE                       !==  BLK formulation  ==!   turbulent fluxes computation
+         
+         !                             ! Compute true air density :
+         IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
+            zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , pslp(:,:) )
+         ELSE                                      ! At zt:
+            zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) )
+         END IF
+         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zztmp         = zrhoa(ji,jj) * zU_zu(ji,jj)
+!!gm to be done by someone: check the efficiency of the call of cp_air within do loops 
+               psen  (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * Ch_atm(ji,jj) * (  zst(ji,jj) - t_zu(ji,jj) )
+               pevp  (ji,jj) = rn_efac*MAX( 0._wp,      zztmp * Ce_atm(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) )
+               zztmp         = zztmp * Cd_atm(ji,jj)
+               pcd_du(ji,jj) = zztmp
+               taum  (ji,jj) = zztmp * wndm  (ji,jj) 
+               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+            END DO
+         END DO
+
+         CALL iom_put( "taum_oce", taum )   ! output wind stress module
+         
+         ! ... utau, vtau at U- and V_points, resp.
+         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1
+               utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
+                  &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
+               vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
+                  &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
+
+         IF(ln_ctl) THEN
+            CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_1: wndm   : ')
+            CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   &
+               &          tab2d_2=vtau  , clinfo2='            vtau   : ', mask2=vmask )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF(ln_ctl) THEN
+         CALL prt_ctl( tab2d_1=pevp  , clinfo1=' blk_oce_1: pevp   : ' )
+         CALL prt_ctl( tab2d_1=psen  , clinfo1=' blk_oce_1: psen   : ' )
+         CALL prt_ctl( tab2d_1=pssq  , clinfo1=' blk_oce_1: pssq   : ' )
+      ENDIF
+      !
+  END SUBROUTINE blk_oce_1
+
+
+   SUBROUTINE blk_oce_2( kt, ptair, pqsr, pqlw, pprec,   &   ! <<= in
+              &              psnow, pst , psen, pevp     )   ! <<= in
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_oce_2  ***
+      !!
+      !! ** Purpose :   finalize the momentum, heat and freshwater fluxes computation 
+      !!                at the ocean surface at each time step knowing Cd, Ch, Ce and 
+      !!                atmospheric variables (from ABL or external data)
       !!
       !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
@@ -360,55 +563,24 @@
       !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
       !!              - emp     : evaporation minus precipitation       (kg/m2/s)
-      !!
-      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
-      !!---------------------------------------------------------------------
-      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
-      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
-      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
-      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
-      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
+      !!---------------------------------------------------------------------
+      INTEGER , INTENT(in)                 ::   kt    ! time step index
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptair
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqsr
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqlw
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pevp
       !
       INTEGER  ::   ji, jj               ! dummy loop indices
-      REAL(wp) ::   zztmp                ! local variable
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst
-      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes
-      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation
+      REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable
+      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes
+      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat fluxes and evaporation
       REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin
-      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
       !!---------------------------------------------------------------------
       !
       ! local scalars ( place there for vector optimisation purposes)
       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
-
-      ! ----------------------------------------------------------------------------- !
-      !      0   Wind components and module at T-point relative to the moving ocean   !
-      ! ----------------------------------------------------------------------------- !
-
-      ! ... components ( U10m - U_oce ) at T-point (unmasked)
-!!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
-      zwnd_i(:,:) = 0._wp
-      zwnd_j(:,:) = 0._wp
-#if defined key_cyclone
-      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
-            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
-         END DO
-      END DO
-#endif
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
-            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
-         END DO
-      END DO
-      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
-      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
-      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
-         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
 
       ! ----------------------------------------------------------------------------- !
@@ -417,103 +589,19 @@
 
       ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
-      zztmp = 1. - albo
-      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
-      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
-      ENDIF
-
-      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
-
-      ! ----------------------------------------------------------------------------- !
-      !     II    Turbulent FLUXES                                                    !
-      ! ----------------------------------------------------------------------------- !
-
-      ! ... specific humidity at SST and IST tmask(
-      zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
-      !!
-      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
-      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
-      !!    (since reanalysis products provide T at z, not theta !)
-      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
-
-      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
-      !
-      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE DEFAULT
-         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
-      END SELECT
-
-      !                          ! Compute true air density :
-      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
-         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
-      ELSE                                      ! At zt:
-         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
-      END IF
-
-!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
-!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
-
-      DO jj = 1, jpj             ! tau module, i and j component
-         DO ji = 1, jpi
-            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
-            taum  (ji,jj) = zztmp * wndm  (ji,jj)
-            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
-            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
-         END DO
-      END DO
-
-      !                          ! add the HF tau contribution to the wind stress module
-      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
-
-      CALL iom_put( "taum_oce", taum )   ! output wind stress module
-
-      ! ... utau, vtau at U- and V_points, resp.
-      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
-      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1
-            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
-               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
-            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
-               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
-         END DO
-      END DO
-      CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
+      zztmp = 1._wp - albo
+      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( pqsr(:,:) ) * tmask(:,:,1)
+      ELSE                  ;   qsr(:,:) = zztmp *          pqsr(:,:)   * tmask(:,:,1)
+      ENDIF
+
+      zqlw(:,:) = ( pqlw(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
 
       !  Turbulent fluxes over ocean
       ! -----------------------------
 
-      ! zqla used as temporary array, for rho*U (common term of bulk formulae):
-      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
-
-      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
-         !! q_air and t_air are given at 10m (wind reference height)
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
-      ELSE
-         !! q_air and t_air are not given at 10m (wind reference height)
-         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
-      ENDIF
-
-      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
-
+      zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:)     ! Latent Heat flux
 
       IF(ln_ctl) THEN
-         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
-         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
-         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
-         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
-         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
-            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
-         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
-         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
+         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_2: zqla   : ' )
+         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_2: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
       ENDIF
 
@@ -522,18 +610,18 @@
       ! ----------------------------------------------------------------------------- !
       !
-      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
-         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
-      !
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
-         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
-         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
-         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
-         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
-         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
-         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
-      qns(:,:) = qns(:,:) * tmask(:,:,1)
+      emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac  ) * tmask(:,:,1)             ! mass flux (evap. - precip.)
+      !
+      zz1 = rn_pfac * rLfus
+      zz2 = rn_pfac * rcp
+      zz3 = rn_pfac * rcpi
+      qns(:,:) = (  zqlw(:,:) - psen(:,:) - zqla(:,:)                          &   ! Downward Non Solar
+         &        - psnow(:,:) * zztmp                                         &   ! remove latent melting heat for solid precip
+         &        - pevp(:,:) * pst(:,:) * rcp                                 &   ! remove evap heat content at SST
+         &        + ( pprec(:,:) - psnow(:,:) ) * ( ptair(:,:) - rt0 ) * zz2   &   ! add liquid precip heat content at Tair
+         &        + psnow(:,:) * ( MIN( ptair(:,:), rt0 ) - rt0 ) * zz3        &   ! add solid  precip heat content at min(Tair,Tsnow)
+         &       ) * tmask(:,:,1)
       !
 #if defined key_si3
-      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3)
+      qns_oce(:,:) = zqlw(:,:) - psen(:,:) - zqla(:,:)                             ! non solar without emp (only needed by SI3)
       qsr_oce(:,:) = qsr(:,:)
 #endif
@@ -541,12 +629,12 @@
       IF ( nn_ice == 0 ) THEN
          CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
-         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
+         CALL iom_put( "qsb_oce" , - psen )                 ! output downward sensible heat over the ocean
          CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
-         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
+         CALL iom_put( "qemp_oce",   qns-zqlw+psen+zqla )   ! output downward heat content of E-P over the ocean
          CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
-         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
-         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
+         tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
+         sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
          CALL iom_put( 'snowpre', sprecip )                 ! Snow
          CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
@@ -554,15 +642,12 @@
       !
       IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
-         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
-         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
-         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
-            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
-      ENDIF
-      !
-   END SUBROUTINE blk_oce
-
-
-
+         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ')
+         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla  : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
+         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ')
+      ENDIF
+      !
+   END SUBROUTINE blk_oce_2
+
+   
    FUNCTION rho_air( ptak, pqa, pslp )
       !!-------------------------------------------------------------------------------
@@ -584,5 +669,22 @@
 
 
-   FUNCTION cp_air( pqa )
+   FUNCTION cp_air_0d( pqa )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air  ***
+      !!
+      !! ** Purpose : Compute specific heat (Cp) of moist air
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!-------------------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
+      REAL(wp)             ::   cp_air_0d! specific heat of moist air   [J/K/kg]
+      !!-------------------------------------------------------------------------------
+      !
+      Cp_air_0d = Cp_dry + Cp_vap * pqa
+      !
+   END FUNCTION cp_air_0d
+
+
+   FUNCTION cp_air_2d( pqa )
       !!-------------------------------------------------------------------------------
       !!                           ***  FUNCTION cp_air  ***
@@ -593,10 +695,10 @@
       !!-------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_2d! specific heat of moist air   [J/K/kg]
       !!-------------------------------------------------------------------------------
       !
-      Cp_air = Cp_dry + Cp_vap * pqa
-      !
-   END FUNCTION cp_air
+      Cp_air_2d = Cp_dry + Cp_vap * pqa
+      !
+   END FUNCTION cp_air_2d
 
 
@@ -714,6 +816,4 @@
       Ce_atm(:,:) = Cd_ice
 
-      wndm_ice(:,:) = 0._wp      !!gm brutal....
-
       ! ------------------------------------------------------------ !
       !    Wind module relative to the moving ice ( U10m - U_ice )   !
@@ -743,9 +843,4 @@
       ! Computing density of air! Way denser that 1.2 over sea-ice !!!
       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
-
-      !!gm brutal....
-      utau_ice  (:,:) = 0._wp
-      vtau_ice  (:,:) = 0._wp
-      !!gm end
 
       ! ------------------------------------------------------------ !
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbccpl.F90	(revision 11273)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbccpl.F90	(revision 11275)
@@ -533,5 +533,4 @@
       !                                                      ! ------------------------- !
       srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
-      lhftau = srcv(jpr_taum)%laction
       !
       !                                                      ! ------------------------- !
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcmod.F90	(revision 11273)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcmod.F90	(revision 11275)
@@ -15,4 +15,5 @@
    !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting
    !!            4.0  ! 2016-06  (L. Brodeau) new general bulk formulation
+   !!            4.0  ! 2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE)
    !!----------------------------------------------------------------------
 
@@ -32,4 +33,5 @@
    USE sbcflx         ! surface boundary condition: flux formulation
    USE sbcblk         ! surface boundary condition: bulk formulation
+   USE sbcabl         ! atmospheric boundary layer
    USE sbcice_if      ! surface boundary condition: ice-if sea-ice model
 #if defined key_si3
@@ -92,10 +94,10 @@
       !!
       NAMELIST/namsbc/ nn_fsbc  ,                                                    &
-         &             ln_usr   , ln_flx   , ln_blk       ,                          &
+         &             ln_usr   , ln_flx   , ln_blk   , ln_abl,                      &
          &             ln_cpl   , ln_mixcpl, nn_components,                          &
          &             nn_ice   , ln_ice_embd,                                       &
          &             ln_traqsr, ln_dm2dc ,                                         &
          &             ln_rnf   , nn_fwb   , ln_ssr   , ln_isf    , ln_apr_dyn ,     &
-         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor   ,     &
+         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor  ,     &
          &             ln_tauw  , nn_lsm, nn_sdrift
       !!----------------------------------------------------------------------
@@ -137,4 +139,5 @@
          WRITE(numout,*) '         flux         formulation                   ln_flx        = ', ln_flx
          WRITE(numout,*) '         bulk         formulation                   ln_blk        = ', ln_blk
+         WRITE(numout,*) '         ABL          formulation                   ln_abl        = ', ln_abl
          WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : '
          WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
@@ -225,7 +228,11 @@
       CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)  
       CASE( 2 )                        !- SI3  ice model
+         IF( .NOT.( ln_blk .OR. ln_cpl .OR. ln_abl ) )   &
+            &                   CALL ctl_stop( 'sbc_init : SI3 sea-ice model requires ln_blk or ln_cpl or ln_abl = T' )
       CASE( 3 )                        !- CICE ice model
-         IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl = T' )
-         IF( lk_agrif                    )   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' ) 
+         IF( .NOT.( ln_blk .OR. ln_cpl .OR. ln_abl ) )   &
+            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl or ln_abl = T' )
+         IF( lk_agrif                                )   &
+            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' ) 
       CASE DEFAULT                     !- not supported
       END SELECT
@@ -270,4 +277,5 @@
       IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation
       IF( ln_blk          ) THEN   ;   nsbc = jp_blk     ; icpt = icpt + 1   ;   ENDIF       ! bulk                 formulation
+      IF( ln_abl          ) THEN   ;   nsbc = jp_abl     ; icpt = icpt + 1   ;   ENDIF       ! ABL                  formulation
       IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
       IF( ll_opa          ) THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
@@ -281,4 +289,5 @@
          CASE( jp_flx     )   ;   WRITE(numout,*) '   ==>>>   flux formulation'
          CASE( jp_blk     )   ;   WRITE(numout,*) '   ==>>>   bulk formulation'
+         CASE( jp_abl     )   ;   WRITE(numout,*) '   ==>>>   ABL  formulation'
          CASE( jp_purecpl )   ;   WRITE(numout,*) '   ==>>>   pure coupled formulation'
 !!gm abusive use of jp_none ??   ===>>> need to be check and changed by adding a jp_sas parameter
@@ -327,4 +336,6 @@
       IF( ln_blk      )   CALL sbc_blk_init            ! bulk formulae initialization
 
+      IF( ln_abl      )    CALL sbc_abl_init           ! Atmospheric Boundary Layer (ABL)
+
       IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization
       !
@@ -355,5 +366,5 @@
          CALL iom_set_rstw_var_active('sfx_b')
       ENDIF
-
+      ! 
    END SUBROUTINE sbc_init
 
@@ -424,4 +435,8 @@
          IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA
                                CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean
+                               !
+      CASE( jp_abl     )
+         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA
+                               CALL sbc_abl       ( kt )                    ! ABL  formulation for the ocean
                                !
       CASE( jp_purecpl )   ;   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! pure coupled formulation
@@ -498,5 +513,5 @@
             CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b, ldxios = lrxios )   ! before i-stress  (U-point)
             CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b, ldxios = lrxios )   ! before j-stress  (V-point)
-            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b, ldxios = lrxios  )   ! before non solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo,  'qns_b',  qns_b, ldxios = lrxios )   ! before non solar heat flux (T-point)
             ! The 3D heat content due to qsr forcing is treated in traqsr
             ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b, ldxios = lrxios  ) ! before     solar heat flux (T-point)
