Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/DOM/phycst.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/DOM/phycst.F90	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/DOM/phycst.F90	(revision 11615)
@@ -36,5 +36,5 @@
    REAL(wp), PUBLIC ::   omega                       !: earth rotation parameter           [s-1]
    REAL(wp), PUBLIC ::   ra       = 6371229._wp      !: earth radius                       [m]
-   REAL(wp), PUBLIC ::   grav     = 9.80665_wp       !: gravity                            [m/s2]   
+   REAL(wp), PARAMETER, PUBLIC ::   grav     = 9.80665_wp       !: gravity                            [m/s2]    !LB
    REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin]
 
@@ -82,4 +82,5 @@
    REAL(wp), PARAMETER, PUBLIC :: rho0_a  = 1.2_wp      !: Approx. of density of air                       [kg/m^3]
    REAL(wp), PARAMETER, PUBLIC :: rho0_w  = 1025._wp    !: Density of sea-water  (ECMWF->1025)             [kg/m^3]
+   REAL(wp), PARAMETER, PUBLIC :: roadrw = rho0_a/rho0_w !: Density ratio
    REAL(wp), PARAMETER, PUBLIC :: rCp0_w  = 4190._wp    !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg]
    REAL(wp), PARAMETER, PUBLIC :: rnu0_w  = 1.e-6_wp    !: kinetic viscosity of water                      [m^2/s]
@@ -90,5 +91,7 @@
    !                                                    !: the small fraction of downwelling shortwave reflected at the
    !                                                    !: surface (Lind & Katsaros, 1986)
-   REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp  !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt  
+   REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp  !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt
+   REAL(wp), PARAMETER, PUBLIC :: rtt0 = 273.16_wp        !: triple point of temperature    [K]
+   REAL(wp), PARAMETER, PUBLIC :: rcst_cs = 16._wp*grav*rho0_w*rCp0_w*rnu0_w*rnu0_w*rnu0_w/(rk0_w*rk0_w) !: for cool-skin parameterizations...
    !#LB.
 
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90	(revision 11615)
@@ -153,5 +153,5 @@
    !!                     Cool-skin/Warm-layer
    !!----------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature (used if ln_skin==.true.)  [K] !LB
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature (used if ln_skin_cs==.true. .OR. ln_skin_wl==.true.)  [K] !LB
 
    
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11615)
@@ -48,5 +48,5 @@
    USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
    USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013)
-   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)
+   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 45r1)
    !
    USE iom            ! I/O manager library
@@ -88,5 +88,5 @@
    LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
    LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013)
-   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
+   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 45r1)
    !
    LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data
@@ -108,5 +108,6 @@
 
    !LB:
-   LOGICAL  ::   ln_skin        ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
+   LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB
+   LOGICAL  ::   ln_skin_wl     ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
    LOGICAL  ::   ln_humi_sph    ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB
    LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
@@ -125,5 +126,5 @@
    INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
    INTEGER, PARAMETER ::   np_COARE_3p6 = 3   ! "COARE 3.6" algorithm   (Edson et al. 2013)
-   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31)
+   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 45r1)
 
    !! * Substitutions
@@ -183,5 +184,5 @@
          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &
          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
-         &                 ln_skin, ln_humi_sph, ln_humi_dpt, ln_humi_rlh    ! cool-skin / warm-layer !LB
+         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh    ! cool-skin / warm-layer !LB
       !!---------------------------------------------------------------------
       !
@@ -221,5 +222,5 @@
       !LB:
       !                             !** initialization of the cool-skin / warm-layer parametrization
-      IF( ln_skin ) THEN
+      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
          IF ( ln_NCAR ) CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' )
          !                       ! allocate array(s) for cool-skin/warm-layer param.
@@ -246,5 +247,5 @@
             CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
                &           '              ==> We force time interpolation = .false. for qsr' )
-            sn_qsr%ln_tint = .FALSE.
+            sn_qsr%ln_tint = .false.
          ENDIF
       ENDIF
@@ -301,5 +302,5 @@
          WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
          WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6
-         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
+         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 45r1)            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
@@ -317,9 +318,10 @@
          CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
          CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
-         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)'
+         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 45r1)'
          END SELECT
          !
          WRITE(numout,*)
-         WRITE(numout,*) '      use cool-skin/warm-layer parameterization (SSST)     ln_skin     = ', ln_skin !LB
+         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs !LB
+         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl !LB
          !
          !LB:
@@ -497,5 +499,5 @@
       zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
 
-      IF ( ln_skin ) THEN
+      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
          zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action
          zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action
@@ -522,28 +524,29 @@
 
 
-      IF ( ln_skin ) THEN
+      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
 
          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
 
          CASE( np_COARE_3p0 )
-            CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
-               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
+            CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                            &  ! COARE v3.0
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
                &                Tsk_b=tsk(:,:) )
 
          CASE( np_COARE_3p6 )
-            CALL turb_coare3p6 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6
-               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
+            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day=', nsec_day
+            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
-               &                Tsk_b=tsk(:,:) )
+               &                isecday_utc=nsec_day, plong=glamt(:,:) )
 
          CASE( np_ECMWF     )
-            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF
-               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
-               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
-               &                Tsk_b=tsk(:,:) )
+            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day=', nsec_day
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   &
+               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
 
          CASE DEFAULT
-            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' )
+            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' )
          END SELECT
 
@@ -562,11 +565,12 @@
          tsk(:,:) = zst(:,:)
 
-      ELSE !IF ( ln_skin )
+      ELSE !IF ( ln_skin_cs .OR. ln_skin_wl )
 
 
          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
             !
-         CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2
-            &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         CASE( np_NCAR      )
+            CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, 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 )
@@ -575,17 +579,19 @@
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
-         CASE( np_COARE_3p6 )   ;   CALL turb_coare3p6( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6
-            &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         CASE( np_COARE_3p6 )
+            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!'
+            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
          CASE( np_ECMWF     )
             IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
-            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! 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
 
-      END IF ! IF ( ln_skin )
+      END IF ! IF ( ln_skin_cs .OR. ln_skin_wl )
 
 
@@ -653,5 +659,5 @@
       ! ----------------------------------------------------------------------------- !
 
-      !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin==.TRUE.)
+      !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
       zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux
 
@@ -702,5 +708,5 @@
       CALL iom_put( 'snowpre', sprecip )                 ! Snow
       CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
-      IF ( ln_skin ) THEN
+      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
          CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
          CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11615)
@@ -2,4 +2,6 @@
    !!======================================================================
    !!                   ***  MODULE  sbcblk_algo_coare3p0  ***
+   !!
+   !!                After Fairall et al, 2003
    !! Computes:
    !!   * bulk transfer coefficients C_D, C_E and C_H
@@ -7,6 +9,4 @@
    !!   * the effective bulk wind speed at 10m U_blk
    !!   => all these are used in bulk formulas in sbcblk.F90
-   !!
-   !!    Using the bulk formulation/param. of COARE v3, Fairall et al., 2003
    !!
    !!       Routine turb_coare3p0 maintained and developed in AeroBulk
@@ -38,10 +38,10 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
-   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+   USE sbcblk_skin_coare3p6  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC ::   TURB_COARE3P0   ! called by sbcblk.F90
+   PUBLIC :: COARE3P0_INIT, TURB_COARE3P0
 
    !                                              !! COARE own values for given constants:
@@ -54,9 +54,46 @@
 CONTAINS
 
-   SUBROUTINE turb_coare3p0( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
-      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
+
+   SUBROUTINE coare3p0_init(l_use_cs, l_use_wl)
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION coare3p0_init  ***
+      !!
+      !! INPUT :
+      !! -------
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
+      !!---------------------------------------------------------------------
+      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization
+      INTEGER :: ierr
+      !!---------------------------------------------------------------------
+      IF ( l_use_wl ) THEN
+         ierr = 0
+         PRINT *, ' *** coare3p0_init: WL => allocating Tau_ac, Qnt_ac, and H_wl :', jpi,jpj
+         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr )
+         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac and Qnt_ac failed!'
+         Tau_ac(:,:) = 0._wp
+         Qnt_ac(:,:)   = 0._wp
+         H_wl(:,:)    = H_wl_max
+         PRINT *, ' *** Tau_ac , Qnt_ac, and H_wl allocated!'
+      END IF
+      !!
+      IF ( l_use_cs ) THEN
+         ierr = 0
+         PRINT *, ' *** coare3p0_init: CS => allocating delta_vl :', jpi,jpj
+         ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr )
+         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of delta_vl and Qnt_ac failed!'
+         delta_vl(:,:) = 0.001_wp      ! First guess of zdelta [m]
+         PRINT *, ' *** delta_vl allocated!'
+      END IF
+   END SUBROUTINE coare3p0_init
+
+
+
+   SUBROUTINE turb_coare3p0( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl,  &
+      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                               &
       &                      Cdn, Chn, Cen,                      &
-      &                      Qsw, rad_lw, slp,                   &
-      &                      Tsk_b )
+      &                      Qsw, rad_lw, slp, pdT_cs,                                    & ! optionals for cool-skin (and warm-layer)
+      &                      isecday_utc, plong, pdT_wl, Hwl )                             ! optionals for warm-layer only
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_coare3p0  ***
@@ -74,9 +111,12 @@
       !! INPUT :
       !! -------
+      !!    *  kt   : current time step (starts at 1)
       !!    *  zt   : height for temperature and spec. hum. of air            [m]
       !!    *  zu   : height for wind speed (usually 10m)                     [m]
-      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
       !!    *  t_zt : potential air temperature at zt                         [K]
       !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
+      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
       !!
       !! INPUT/OUTPUT:
@@ -87,13 +127,17 @@
       !!
       !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
-      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True)
-      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False)
-      !!
-      !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!):
+      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True)
+      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False)
+      !!
+      !! OPTIONAL INPUT:
       !! ---------------
       !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
       !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
       !!    *  slp    : sea-level pressure                                    [Pa]
-      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
+      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
+      !!    * isecday_utc:
+      !!    *  plong  : longitude array                                       [deg.E]
+      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
+      !!    * Hwl     : depth of warm layer                                   [m]
       !!
       !! OUTPUT :
@@ -109,4 +153,5 @@
       !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
+      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
       REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
@@ -116,4 +161,6 @@
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
+      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
@@ -127,5 +174,10 @@
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
-      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
+      !
+      INTEGER,  INTENT(in   ), OPTIONAL                     ::   isecday_utc ! current UTC time, counted in second since 00h of the current day
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   plong    !             [deg.E]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m]
       !
       INTEGER :: j_itt
@@ -141,31 +193,49 @@
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
       !
-      ! Cool skin:
-      LOGICAL :: l_use_skin = .FALSE.
-      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
+         &                zsst,   &  ! to back up the initial bulk SST
+         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
+         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K]
+         &                zHwl       ! depth of warm-layer [m]
+
+      !
+      LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE.
+      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0.F90'
+      CHARACTER(len=128) :: cf_tmp
       !!----------------------------------------------------------------------------------
 
-      ! Cool skin ?
-      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
-         l_use_skin = .TRUE.
-      END IF
-      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
+      IF ( kt == 1 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays
+
 
       l_zt_equal_zu = .FALSE.
       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
-
       IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
 
-      !! Initialization for cool skin:
-      zsst   = T_s    ! save the bulk SST
-      IF( l_use_skin ) THEN
-         ! First guess for skin temperature:
-         IF( PRESENT(Tsk_b) ) THEN
-            T_s = Tsk_b
-         ELSE
-            T_s = T_s - 0.25     ! sst - 0.25
-         END IF
-         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
-      END IF
+      !! Initializations for cool skin and warm layer:
+      IF ( l_use_cs ) THEN
+         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN
+            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTc(jpi,jpj) )
+         pdTc(:,:) = -0.25_wp  ! First guess of skin correction
+      END IF
+
+      IF ( l_use_wl ) THEN
+         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN
+            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTw(jpi,jpj) )
+         IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) )
+      END IF
+
+      IF ( l_use_cs .OR. l_use_wl ) THEN
+         ALLOCATE ( zsst(jpi,jpj) )
+         zsst = T_s ! backing up the bulk SST
+         IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!!
+      END IF
+
 
       !! First guess of temperature and humidity at height zu:
@@ -190,6 +260,5 @@
       z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
-      ztmp2  = vkarmn/ztmp0
-      Cd     = ztmp2*ztmp2    ! first guess of Cd
+      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
 
       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
@@ -234,9 +303,8 @@
          ztmp1 = u_star*u_star   ! u*^2
 
-         !! Update wind at zu taking into acount convection-related wind gustiness:
-         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
-         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2
-         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
-         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
+         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8):
+         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2
+         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
+         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
 
@@ -251,10 +319,10 @@
          !! Adjustment the wind at 10m (not needed in the current algo form):
          !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
-         
+
          !! Roughness lengthes z0, z0t (z0q = z0t) :
          ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
          z0    = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6)
-         ztmp1 = z0*u_star/znu_a                                           ! Re_r: roughness Reynolds number
-         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!!
+         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp     ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
+         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!!
 
          !! Turbulent scales at zu :
@@ -267,31 +335,51 @@
 
          IF( .NOT. l_zt_equal_zu ) THEN
-            ! What's need to be done if zt /= zu
-            !! Re-updating temperature and humidity at zu :
-            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
-            ztmp1 = log(zt/zu) + ztmp2
+            !! Re-updating temperature and humidity at zu if zt /= zu :
+            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t)
             t_zu = t_zt - t_star/vkarmn*ztmp1
             q_zu = q_zt - q_star/vkarmn*ztmp1
          END IF
 
-         !! SKIN related part
-         !! -----------------
-         IF( l_use_skin ) THEN
-            !! compute transfer coefficients at zu : lolo: verifier...
-            ztmp0 = u_star/U_blk
-            Ch   = ztmp0*t_star/dt_zu
-            Ce   = ztmp0*q_star/dq_zu
-            ! Non-Solar heat flux to the ocean:
-            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
-            ztmp2 = T_s*T_s
-            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
-               &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux
-            !!         => "ztmp1" is the net non-solar surface heat flux !
-            !! Updating the values of the skin temperature T_s and q_s :
-            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
-            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given
-         END IF
-
-         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
+
+         IF( l_use_cs ) THEN
+            !! Cool-skin contribution
+
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
+
+            CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
+
+            T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1)
+            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+
+         END IF
+
+         IF( l_use_wl ) THEN
+            !! Warm-layer contribution
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
+            !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
+            IF (PRESENT(Hwl)) THEN
+               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
+            ELSE
+               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw )
+            END IF
+            !! Updating T_s and q_s !!!
+            T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1)
+            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+            !LOLO:
+            !IF ( j_itt == nb_itt) THEN
+            !   WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea
+            !   CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac')
+            !   WRITE(cf_tmp,  '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea
+            !   CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw')
+            !END IF
+            !LOLO.
+         END IF
+
+
+         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
             dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
             dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
@@ -312,4 +400,14 @@
 
       IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
+
+      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1)
+      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1)
+
+      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+      IF (          l_use_cs      ) DEALLOCATE ( pdTc )
+      IF (          l_use_wl      ) THEN
+         DEALLOCATE ( pdTw )
+         IF (PRESENT(Hwl)) DEALLOCATE ( zHwl )
+      END IF
 
    END SUBROUTINE turb_coare3p0
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11615)
@@ -38,10 +38,10 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
-   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+   USE sbcblk_skin_coare3p6  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC ::   TURB_COARE3P6   ! called by sbcblk.F90
+   PUBLIC :: COARE3P6_INIT, TURB_COARE3P6
 
    !                                              !! COARE own values for given constants:
@@ -54,9 +54,46 @@
 CONTAINS
 
-   SUBROUTINE turb_coare3p6( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
-      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
+
+   SUBROUTINE coare3p6_init(l_use_cs, l_use_wl)
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION coare3p6_init  ***
+      !!
+      !! INPUT :
+      !! -------
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
+      !!---------------------------------------------------------------------
+      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization
+      INTEGER :: ierr
+      !!---------------------------------------------------------------------
+      IF ( l_use_wl ) THEN
+         ierr = 0
+         PRINT *, ' *** coare3p6_init: WL => allocating Tau_ac, Qnt_ac, and H_wl :', jpi,jpj
+         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr )
+         !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of Tau_ac and Qnt_ac failed!'
+         Tau_ac(:,:) = 0._wp
+         Qnt_ac(:,:)   = 0._wp
+         H_wl(:,:)    = H_wl_max
+         PRINT *, ' *** Tau_ac , Qnt_ac, and H_wl allocated!'
+      END IF
+      !!
+      IF ( l_use_cs ) THEN
+         ierr = 0
+         PRINT *, ' *** coare3p6_init: CS => allocating delta_vl :', jpi,jpj
+         ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr )
+         !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of delta_vl and Qnt_ac failed!'
+         delta_vl(:,:) = 0.001_wp      ! First guess of zdelta [m]
+         PRINT *, ' *** delta_vl allocated!'
+      END IF
+   END SUBROUTINE coare3p6_init
+
+
+
+   SUBROUTINE turb_coare3p6( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl,  &
+      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                               &
       &                      Cdn, Chn, Cen,                      &
-      &                      Qsw, rad_lw, slp,                   &
-      &                      Tsk_b )
+      &                      Qsw, rad_lw, slp, pdT_cs,                                    & ! optionals for cool-skin (and warm-layer)
+      &                      isecday_utc, plong, pdT_wl, Hwl )                             ! optionals for warm-layer only
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_coare3p6  ***
@@ -74,9 +111,12 @@
       !! INPUT :
       !! -------
+      !!    *  kt   : current time step (starts at 1)
       !!    *  zt   : height for temperature and spec. hum. of air            [m]
       !!    *  zu   : height for wind speed (usually 10m)                     [m]
-      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
       !!    *  t_zt : potential air temperature at zt                         [K]
       !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
+      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
       !!
       !! INPUT/OUTPUT:
@@ -90,10 +130,14 @@
       !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False)
       !!
-      !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!):
+      !! OPTIONAL INPUT:
       !! ---------------
       !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
       !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
       !!    *  slp    : sea-level pressure                                    [Pa]
-      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
+      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
+      !!    * isecday_utc:
+      !!    *  plong  : longitude array                                       [deg.E]
+      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
+      !!    * Hwl     : depth of warm layer                                   [m]
       !!
       !! OUTPUT :
@@ -109,4 +153,5 @@
       !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
+      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
       REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
@@ -116,4 +161,6 @@
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
+      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
@@ -127,5 +174,10 @@
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
-      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
+      !
+      INTEGER,  INTENT(in   ), OPTIONAL                     ::   isecday_utc ! current UTC time, counted in second since 00h of the current day
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   plong    !             [deg.E]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m]
       !
       INTEGER :: j_itt
@@ -141,31 +193,49 @@
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
       !
-      ! Cool skin:
-      LOGICAL :: l_use_skin = .FALSE.
-      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
+         &                zsst,   &  ! to back up the initial bulk SST
+         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
+         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K]
+         &                zHwl       ! depth of warm-layer [m]
+
+      !
+      LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE.
+      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p6.F90'
+      CHARACTER(len=128) :: cf_tmp
       !!----------------------------------------------------------------------------------
 
-      ! Cool skin ?
-      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
-         l_use_skin = .TRUE.
-      END IF
-      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
+      IF ( kt == 1 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays
+
 
       l_zt_equal_zu = .FALSE.
       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
-
       IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
 
-      !! Initialization for cool skin:
-      zsst   = T_s    ! save the bulk SST
-      IF( l_use_skin ) THEN
-         ! First guess for skin temperature:
-         IF( PRESENT(Tsk_b) ) THEN
-            T_s = Tsk_b
-         ELSE
-            T_s = T_s - 0.25     ! sst - 0.25
-         END IF
-         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
-      END IF
+      !! Initializations for cool skin and warm layer:
+      IF ( l_use_cs ) THEN
+         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN
+            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTc(jpi,jpj) )
+         pdTc(:,:) = -0.25_wp  ! First guess of skin correction
+      END IF
+
+      IF ( l_use_wl ) THEN
+         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN
+            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTw(jpi,jpj) )
+         IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) )
+      END IF
+
+      IF ( l_use_cs .OR. l_use_wl ) THEN
+         ALLOCATE ( zsst(jpi,jpj) )
+         zsst = T_s ! backing up the bulk SST
+         IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!!
+      END IF
+
 
       !! First guess of temperature and humidity at height zu:
@@ -190,6 +260,5 @@
       z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
-      ztmp2  = vkarmn/ztmp0
-      Cd     = ztmp2*ztmp2    ! first guess of Cd
+      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
 
       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
@@ -234,9 +303,8 @@
          ztmp1 = u_star*u_star   ! u*^2
 
-         !! Update wind at zu taking into acount convection-related wind gustiness:
-         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
-         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2
-         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
-         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
+         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8):
+         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2
+         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
+         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
 
@@ -251,10 +319,10 @@
          !! Adjustment the wind at 10m (not needed in the current algo form):
          !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
-         
+
          !! Roughness lengthes z0, z0t (z0q = z0t) :
          ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
-         z0    = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6)
-         ztmp1 = z0*u_star/znu_a                                           ! Re_r: roughness Reynolds number
-         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1**(-0.72_wp))   ! COARE 3.6
+         z0    = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ]
+         ztmp1 = ( znu_a / (z0*u_star) )**0.72_wp     ! COARE3.6-specific! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
+         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1 )   ! COARE3.6-specific!
 
          !! Turbulent scales at zu :
@@ -267,31 +335,51 @@
 
          IF( .NOT. l_zt_equal_zu ) THEN
-            ! What's need to be done if zt /= zu
-            !! Re-updating temperature and humidity at zu :
-            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
-            ztmp1 = log(zt/zu) + ztmp2
+            !! Re-updating temperature and humidity at zu if zt /= zu :
+            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t)
             t_zu = t_zt - t_star/vkarmn*ztmp1
             q_zu = q_zt - q_star/vkarmn*ztmp1
          END IF
 
-         !! SKIN related part
-         !! -----------------
-         IF( l_use_skin ) THEN
-            !! compute transfer coefficients at zu : lolo: verifier...
-            ztmp0 = u_star/U_blk
-            Ch   = ztmp0*t_star/dt_zu
-            Ce   = ztmp0*q_star/dq_zu
-            ! Non-Solar heat flux to the ocean:
-            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
-            ztmp2 = T_s*T_s
-            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
-               &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux
-            !!         => "ztmp1" is the net non-solar surface heat flux !
-            !! Updating the values of the skin temperature T_s and q_s :
-            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
-            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given
-         END IF
-
-         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
+
+         IF( l_use_cs ) THEN
+            !! Cool-skin contribution
+
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
+
+            CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
+
+            T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1)
+            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+
+         END IF
+
+         IF( l_use_wl ) THEN
+            !! Warm-layer contribution
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
+            !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
+            IF (PRESENT(Hwl)) THEN
+               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
+            ELSE
+               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw )
+            END IF
+            !! Updating T_s and q_s !!!
+            T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1)
+            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+            !LOLO:
+            !IF ( j_itt == nb_itt) THEN
+            !   WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea
+            !   CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac')
+            !   WRITE(cf_tmp,  '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea
+            !   CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw')
+            !END IF
+            !LOLO.
+         END IF
+
+
+         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
             dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
             dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
@@ -312,4 +400,14 @@
 
       IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
+
+      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1)
+      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1)
+
+      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+      IF (          l_use_cs      ) DEALLOCATE ( pdTc )
+      IF (          l_use_wl      ) THEN
+         DEALLOCATE ( pdTw )
+         IF (PRESENT(Hwl)) DEALLOCATE ( zHwl )
+      END IF
 
    END SUBROUTINE turb_coare3p6
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11615)
@@ -39,5 +39,5 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
-   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+   USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB
 
    IMPLICIT NONE
@@ -60,14 +60,14 @@
 CONTAINS
 
-   SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
+   SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
       &                      Cdn, Chn, Cen,                      &
-      &                      Qsw, rad_lw, slp,                   &
-      &                      Tsk_b )
+      &                      Qsw, rad_lw, slp, pdT_cs,           & ! optionals for cool-skin (and warm-layer)
+      &                      pdT_wl )                              ! optionals for warm-layer only
       !!----------------------------------------------------------------------
-       !!                      ***  ROUTINE  turb_ecmwf  ***
+      !!                      ***  ROUTINE  turb_ecmwf  ***
       !!
       !! ** Purpose :   Computes turbulent transfert coefficients of surface
-      !!                fluxes according to IFS doc. (cycle 40)
+      !!                fluxes according to IFS doc. (cycle 45r1)
       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
       !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
@@ -82,7 +82,9 @@
       !!    *  zt   : height for temperature and spec. hum. of air            [m]
       !!    *  zu   : height for wind speed (usually 10m)                     [m]
-      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
       !!    *  t_zt : potential air temperature at zt                         [K]
       !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
+      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
       !!
       !! INPUT/OUTPUT:
@@ -101,5 +103,6 @@
       !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
       !!    *  slp    : sea-level pressure                                    [Pa]
-      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
+      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
+      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
       !!
       !! OUTPUT :
@@ -122,4 +125,6 @@
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
+      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
@@ -133,5 +138,7 @@
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
-      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
+      !
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
       !
       INTEGER :: j_itt
@@ -145,35 +152,46 @@
          &  z0, z0t, z0q
       !
-      ! Cool skin:
-      LOGICAL :: l_use_skin = .FALSE.
-      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
+         &                zsst,   &  ! to back up the initial bulk SST
+         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
+         &                pdTw       ! SST increment "dT" for warm layer correction          [K]
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h
       REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
-      !!----------------------------------------------------------------------------------
-
-      ! Cool skin ?
-      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
-         l_use_skin = .TRUE.
-      END IF
-      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
-
-      ! Identical first gess as in COARE, with IFS parameter values though...
-      !
+      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
+      !!----------------------------------------------------------------------------------
+
       l_zt_equal_zu = .FALSE.
       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
 
-      !! Initialization for cool skin:
-      zsst   = T_s    ! save the bulk SST
-      IF( l_use_skin ) THEN
-         ! First guess for skin temperature:
-         IF( PRESENT(Tsk_b) ) THEN
-            T_s = Tsk_b
-         ELSE
-            T_s = T_s - 0.25     ! sst - 0.25
-         END IF
-         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
+      !! Initializations for cool skin and warm layer:
+      IF ( l_use_cs ) THEN
+         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN
+            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTc(jpi,jpj) )
+         pdTc(:,:) = -0.25_wp  ! First guess of skin correction
       END IF
 
+      IF ( l_use_wl ) THEN
+         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN
+            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTw(jpi,jpj) )
+         pdTw(:,:) = 0._wp
+      END IF
+
+      IF ( l_use_cs .OR. l_use_wl ) THEN
+         ALLOCATE ( zsst(jpi,jpj) )
+         zsst = T_s ! backing up the bulk SST
+         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!!
+      END IF
+
+
+      ! Identical first gess as in COARE, with IFS parameter values though...
+      !
       !! First guess of temperature and humidity at height zu:
       t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
@@ -197,6 +215,5 @@
       z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
-      ztmp2  = vkarmn/ztmp0
-      Cd     = ztmp2*ztmp2    ! first guess of Cd
+      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
 
       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
@@ -265,10 +282,8 @@
          z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
 
-         !! Update wind at 10m taking into acount convection-related wind gustiness:
-         !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8
-         ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!!
-         ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn , 0._wp))**(2._wp/3._wp) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
-         !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..)
-         U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp )              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1
+         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8)
+         ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
+         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
+         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
 
@@ -303,23 +318,32 @@
          func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
 
-         !! SKIN related part
-         !! -----------------
-         IF( l_use_skin ) THEN
-            !! compute transfer coefficients at zu : lolo: verifier...
-            Ch = vkarmn*vkarmn/(func_m*func_h)
-            ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q
-            Ce = vkarmn*vkarmn/(func_m*ztmp1)
-            ! Non-Solar heat flux to the ocean:
-            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
-            ztmp2 = T_s*T_s
-            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
-               &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux
-            !!         => "ztmp1" is the net non-solar surface heat flux !
-            !! Updating the values of the skin temperature T_s and q_s :
-            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s )
-            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given
-         END IF
-
-         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
+
+         IF( l_use_cs ) THEN
+            !! Cool-skin contribution
+
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0
+
+            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst, pdTc )  ! Qnsol -> ztmp1
+
+            T_s(:,:) = zsst(:,:) + pdTc(:,:)
+            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+
+         END IF
+
+         IF( l_use_wl ) THEN
+            !! Warm-layer contribution
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2
+            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst, pdTw )
+            !! Updating T_s and q_s !!!
+            T_s(:,:) = zsst(:,:) + pdTw(:,:)
+            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+         END IF
+
+
+         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
             dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
             dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
@@ -337,5 +361,12 @@
       Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
 
-   END SUBROUTINE TURB_ECMWF
+      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc
+      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw
+
+      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+      IF (          l_use_cs      ) DEALLOCATE ( pdTc )
+      IF (          l_use_wl      ) DEALLOCATE ( pdTw )
+
+   END SUBROUTINE turb_ecmwf
 
 
@@ -433,7 +464,4 @@
 
 
-
-
-
    !!======================================================================
 END MODULE sbcblk_algo_ecmwf
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11615)
@@ -29,4 +29,26 @@
    END INTERFACE gamma_moist
 
+   INTERFACE e_sat
+      MODULE PROCEDURE e_sat_vctr, e_sat_sclr
+   END INTERFACE e_sat
+
+   INTERFACE L_vap
+      MODULE PROCEDURE L_vap_vctr, L_vap_sclr
+   END INTERFACE L_vap
+
+   INTERFACE rho_air
+      MODULE PROCEDURE rho_air_vctr, rho_air_sclr
+   END INTERFACE rho_air
+
+   INTERFACE cp_air
+      MODULE PROCEDURE cp_air_vctr, cp_air_sclr
+   END INTERFACE cp_air
+
+      INTERFACE alpha_sw
+      MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr
+   END INTERFACE alpha_sw
+
+
+
    PUBLIC virt_temp
    PUBLIC rho_air
@@ -39,4 +61,7 @@
    PUBLIC q_sat
    PUBLIC q_air_rh
+   !:
+   PUBLIC update_qnsol_tau
+   PUBLIC alpha_sw
 
    !!----------------------------------------------------------------------
@@ -72,7 +97,7 @@
    END FUNCTION virt_temp
 
-   FUNCTION rho_air( ptak, pqa, pslp )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION rho_air  ***
+   FUNCTION rho_air_vctr( ptak, pqa, pslp )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION rho_air_vctr  ***
       !!
       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
@@ -83,10 +108,26 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
-      !!-------------------------------------------------------------------------------
-      !
-      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
-      !
-   END FUNCTION rho_air
+      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air_vctr   ! density of moist air   [kg/m^3]
+      !!-------------------------------------------------------------------------------
+      rho_air_vctr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
+   END FUNCTION rho_air_vctr
+
+   FUNCTION rho_air_sclr( ptak, pqa, pslp )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION rho_air_sclr  ***
+      !!
+      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
+      !!
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!-------------------------------------------------------------------------------
+      REAL(wp), INTENT(in) :: ptak           ! air temperature             [K]
+      REAL(wp), INTENT(in) :: pqa            ! air specific humidity   [kg/kg]
+      REAL(wp), INTENT(in) :: pslp           ! pressure in                [Pa]
+      REAL(wp)             :: rho_air_sclr   ! density of moist air   [kg/m^3]
+      !!-------------------------------------------------------------------------------
+      rho_air_sclr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
+   END FUNCTION rho_air_sclr
+   
+
 
    FUNCTION visc_air(ptak)
@@ -113,23 +154,40 @@
    END FUNCTION visc_air
 
-   FUNCTION L_vap( psst )
+   FUNCTION L_vap_vctr( psst )
       !!---------------------------------------------------------------------------------
-      !!                           ***  FUNCTION L_vap  ***
-      !!
-      !! ** Purpose : Compute the latent heat of vaporization of water out of temperature
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
+      !!                           ***  FUNCTION L_vap_vctr  ***
+      !!
+      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap_vctr   ! latent heat of vaporization   [J/kg]
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
       !!----------------------------------------------------------------------------------
       !
-      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
-      !
-   END FUNCTION L_vap
-
-   FUNCTION cp_air( pqa )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION cp_air  ***
+      L_vap_vctr = (  2.501_wp - 0.00237_wp * ( psst(:,:) - rt0)  ) * 1.e6_wp
+      !
+   END FUNCTION L_vap_vctr
+
+   FUNCTION L_vap_sclr( psst )
+      !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION L_vap_sclr  ***
+      !!
+      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp)             ::   L_vap_sclr   ! latent heat of vaporization   [J/kg]
+      REAL(wp), INTENT(in) ::   psst         ! water temperature                [K]
+      !!----------------------------------------------------------------------------------
+      !
+      L_vap_sclr = (  2.501_wp - 0.00237_wp * ( psst - rt0)  ) * 1.e6_wp
+      !
+   END FUNCTION L_vap_sclr
+
+
+   FUNCTION cp_air_vctr( pqa )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air_vctr  ***
       !!
       !! ** Purpose : Compute specific heat (Cp) of moist air
@@ -138,10 +196,26 @@
       !!-------------------------------------------------------------------------------
       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]
-      !!-------------------------------------------------------------------------------
-      !
-      cp_air = rCp_dry + rCp_vap * pqa
-      !
-   END FUNCTION cp_air
+      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_vctr   ! specific heat of moist air   [J/K/kg]
+      !!-------------------------------------------------------------------------------
+      cp_air_vctr = rCp_dry + rCp_vap * pqa
+   END FUNCTION cp_air_vctr
+
+   FUNCTION cp_air_sclr( pqa )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air_sclr  ***
+      !!
+      !! ** Purpose : Compute specific heat (Cp) of moist air
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!-------------------------------------------------------------------------------
+      REAL(wp), INTENT(in) :: pqa           ! air specific humidity         [kg/kg]
+      REAL(wp)             :: cp_air_sclr   ! specific heat of moist air   [J/K/kg]
+      !!-------------------------------------------------------------------------------
+      cp_air_sclr = rCp_dry + rCp_vap * pqa
+   END FUNCTION cp_air_sclr
+
+
+
+
 
    FUNCTION gamma_moist_vctr( ptak, pqa )
@@ -275,5 +349,37 @@
 
 
-   FUNCTION e_sat_sclr( ptak, pslp )
+   FUNCTION e_sat_vctr(ptak)
+      !!**************************************************
+      !! ptak:     air temperature [K]
+      !! e_sat:  water vapor at saturation [Pa]
+      !!
+      !! Recommended by WMO
+      !!
+      !! Goff, J. A., 1957: Saturation pressure of water on the new kelvin
+      !! temperature scale. Transactions of the American society of heating
+      !! and ventilating engineers, 347–354.
+      !!
+      !! rt0 should be 273.16 (triple point of water) and not 273.15 like here
+      !!**************************************************
+
+      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak    !: temperature (K)
+
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp
+
+      ALLOCATE ( ztmp(jpi,jpj) )
+
+      ztmp(:,:) = rtt0/ptak(:,:)
+
+      e_sat_vctr = 100.*( 10.**(10.79574*(1. - ztmp) - 5.028*LOG10(ptak/rtt0)         &
+         &       + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak/rtt0 - 1.)) )   &
+         &       + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
+
+      DEALLOCATE ( ztmp )
+
+   END FUNCTION e_sat_vctr
+
+
+   FUNCTION e_sat_sclr( ptak )
       !!----------------------------------------------------------------------------------
       !!                   ***  FUNCTION e_sat_sclr  ***
@@ -287,5 +393,4 @@
       !!----------------------------------------------------------------------------------
       REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K]
-      REAL(wp), INTENT(in) ::   pslp    ! sea level atmospheric pressure   [Pa]
       REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg]
       !
@@ -325,5 +430,5 @@
          DO ji = 1, jpi
             !
-            ze_sat =  e_sat_sclr( ptak(ji,jj), pslp(ji,jj) )
+            ze_sat =  e_sat_sclr( ptak(ji,jj) )
             !
             q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat )
@@ -351,5 +456,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj))
+            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
             q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze)
          END DO
@@ -357,4 +462,91 @@
       !
    END FUNCTION q_air_rh
+
+
+   SUBROUTINE UPDATE_QNSOL_TAU( pTs, pqs, pTa, pqa, pust, ptst, pqst, pUb, pslp, prlw, &
+      &                         pQns, pTau,  &
+      &                         Qlat)
+      !!----------------------------------------------------------------------------------
+      !! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw"
+      !!          and the module of the wind stress => pTau = Tau
+      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! air temperature at z=zu [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=zu [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u*
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t*
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqst ! q*
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=zu [m/s]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2]
+      !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
+      !
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat
+      !
+      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zUrho, zTs2, zz0, &
+         &        zQlat, zQsen, zQlw
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      !!----------------------------------------------------------------------------------     
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
+            zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
+            zz0 = pust(ji,jj)/pUb(ji,jj)
+            zCd = zz0*zz0
+            zCh = zz0*ptst(ji,jj)/zdt
+            zCe = zz0*pqst(ji,jj)/zdq
+
+            zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp)     ! rho*U10
+            zTs2  = pTs(ji,jj)*pTs(ji,jj)
+
+            ! Wind stress module:
+            pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo?
+
+            ! Non-Solar heat flux to the ocean:
+            zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp )  ! we do not want a Qlat > 0 !
+            zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt
+            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
+
+            pQns(ji,jj) = zQlat + zQsen + zQlw
+            
+            IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat            
+         END DO
+      END DO
+   END SUBROUTINE UPDATE_QNSOL_TAU
+
+
+   FUNCTION alpha_sw_vctr( psst )
+      !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION alpha_sw_vctr  ***
+      !!
+      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! latent heat of vaporization   [J/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
+      !!----------------------------------------------------------------------------------
+      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
+   END FUNCTION alpha_sw_vctr
+
+   FUNCTION alpha_sw_sclr( psst )
+      !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION alpha_sw_sclr  ***
+      !!
+      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp)             ::   alpha_sw_sclr   ! latent heat of vaporization   [J/kg]
+      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K]
+      !!----------------------------------------------------------------------------------
+      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
+   END FUNCTION alpha_sw_sclr
+
 
 
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin.F90	(revision 11614)
+++ 	(revision )
@@ -1,223 +1,0 @@
-MODULE sbcblk_skin
-   !!======================================================================
-   !!                   ***  MODULE  sbcblk_skin  ***
-   !! Computes:
-   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
-   !!     scheme used at ECMWF
-   !!    Using formulation/param. of IFS of ECMWF (cycle 40r1)
-   !!   based on IFS doc (avaible online on the ECMWF's website)
-   !!
-   !!   From routine turb_ecmwf maintained and developed in AeroBulk
-   !!            (https://github.com/brodeau/aerobulk)
-   !!
-   !! ** Author: L. Brodeau, june 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
-   !!----------------------------------------------------------------------
-   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
-   !!               based on the cool-skin/warm-layer scheme used at ECMWF
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-
-   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
-
-   USE lib_mpp         ! distribued memory computing library
-   USE in_out_manager  ! I/O manager
-   USE lib_fortran     ! to use key_nosignedzero
-
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC :: CSWL_ECMWF   ! called by sbcblk_algo_coare.F90 and sbcblk_algo_ecmwf.F90
-
-   !! Cool-Skin / Warm-Layer related parameters:
-   !LB: we use "rdt" instead !!! REAL(wp), PARAMETER :: rdt0 = 3600.*1.5 !: time step
-   REAL(wp), PARAMETER :: rd0  = 3.        !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005)
-   REAL(wp), PARAMETER :: rNu0 = 0.5       !: Nu (exponent of temperature profile) Eq.11
-   !                                       !: (Zeng & Beljaars 2005) !: set to 0.5 instead of
-   !                                       !: 0.3 to respect a warming of +3 K in calm
-   !                                       !: condition for the insolation peak of +1000W/m^2
-   
-   INTEGER, PARAMETER :: nb_itt_wl = 10    !: number of sub-itterations for solving the differential equation in warm-layer part
-   !                                       !:  => use "nb_itt_wl = 1" for No itteration! => way cheaper !!!
-   !                                       !:    => assumes balance between the last 2 terms of Eq.11 (lhs of eq.11 = 0)
-   !                                       !:    => in that case no need for sub-itterations !
-   !                                       !:    => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions)
-   !                                       !:  => Otherwize use "nb_itt_wl = 10"
-   !
-   !!----------------------------------------------------------------------
-CONTAINS
-
-
-   SUBROUTINE CSWL_ECMWF( pQsw, pQnsol, pustar, pSST, pTs )
-      !!---------------------------------------------------------------------
-      !!
-      !!  Cool-Skin Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
-      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
-      !!
-      !!    As included in IFS Cy40   /  E.C.M.W.F.
-      !!     ------------------------------------------------------------------
-      !!
-      !!  **   INPUT:
-      !!
-      !!     *pQsw*       net solar radiative flux to the ocean
-      !!     *pQnsol*     net non-solar heat flux to the ocean
-      !!     *pustar*     friction velocity u*
-      !!     *pSST*       SST
-      !!
-      !!   **  INPUT/OUTPUT:
-      !!     *pTs*  : as input  =>  previous estimate of skin temperature
-      !!             as output =>  new estimate of skin temperature
-      !!
-      !!------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)    :: pQsw
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)    :: pQnsol
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)    :: pustar
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)    :: pSST
-      !
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: pTs
-      !
-      INTEGER :: ji,jj
-      !
-      REAL(wp) :: rmult, &
-         & zRhoCp_w, &
-         & ZCON2,ZCON3,ZCON4,ZCON5, zQnsol ,zQnet, zlamb, zdelta,&
-         & ZSRD,ZDSST,ZZ,ZEPDU2,&
-         & ZFI,zdL,zdL2, ztmp, &
-         & ZPHI,ZROADRW, &
-         & zus_a, &
-         & zdt, zfs, zsgn
-      !
-      REAL(wp), DIMENSION(jpi,jpj) :: &
-         &  zalpha_w, &       !: thermal expansion coefficient of seawater
-         & zus_w, zus_w2, &   !: u* and u*^2 in water
-         &       zdT_c,   &   !: cool skin temperature increment !lolo rm!
-         &       zdT_w        !: warm skin temperature increment
-      !
-      INTEGER :: nbi, jwl
-      !!------------------------------------------------------------------
-      !
-      !     1. Initialize constants for ocean warm layer and cool skin
-      !
-      !     1.1 General
-      !
-      ZEPDU2  = 0.01_wp   !    security constant for velocity**2   (m2/s2)
-      ZROADRW = rho0_a/rho0_w          ! Density ratio                      (-)
-      zRhoCp_w = rho0_w*rCp0_w
-      !
-      !     1.2C Warm layer parametrization constants
-      !
-      !    ZFI = Fraction of solar radiation absorbed in warm layer (-)
-      ZFI = 1._wp -0.28_wp*EXP(-71.5_wp*rd0) -0.27_wp*EXP(-2.8_wp*rd0) - 0.45_wp*EXP(-0.07_wp*rd0)  !: Eq. 8.135
-      !
-      ZCON3 = rd0*vkarmn*grav/(ZROADRW)**1.5_wp
-      ZCON4 = (rNu0 + 1._wp)*vkarmn/rd0
-      ZCON5 = (rNu0 + 1._wp)/(rNu0*rd0)
-      !
-      !     1.3 Cool skin parametrization constants
-      ZCON2 = 16._wp*grav*zRhoCp_w*rnu0_w**3/(rk0_w*rk0_w)
-      !
-      ! Friction velocities
-      ! "MAX( pustar(:,:), 1.E-4)" is u* in the air !
-      zus_w(:,:)  = MAX( pustar(:,:), 1.E-4_wp)*SQRT(ZROADRW)       ! u* in the water
-      zus_w2(:,:) = zus_w(:,:)*zus_w(:,:)
-      !
-      ! Ocean buoyancy
-      zalpha_w(:,:) = MAX( 1.E-5_wp , 1.E-5_wp*(pTs(:,:) - rt0) ) ! thermal expansion coefficient of water
-      !
-      zdT_c = 0._wp
-      zdT_w = 0._wp
-      !
-      !  3. Cool skin (Fairall et al. 1996)
-      !------------------------------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            ! Non-solar heat loss to the atmosphere:
-            zQnsol = MAX( 1._wp , - pQnsol(ji,jj) )
-
-            zlamb = 6._wp*(1._wp + (zQnsol*zalpha_w(ji,jj)*ZCON2/(zus_w2(ji,jj)*zus_w2(ji,jj)))**0.75)**(-1._wp/3._wp)
-
-            zdelta = zlamb*rnu0_w/zus_w(ji,jj)
-
-            !    Solar absorption
-            zfs   = 0.065_wp + 11._wp*zdelta - (6.6E-5_wp/zdelta)*(1._wp - EXP(-zdelta/8.E-4_wp))  ! Eq. 8.131 / IFS cy40r1, doc, Part IV,
-            zfs   = MAX(zfs , 0.01_wp)
-            zQnet = MAX( 1._wp , -zfs*pQsw(ji,jj) + zQnsol )
-            zdT_c(ji,jj) = -zdelta*zQnet/rk0_w
-
-         END DO
-      END DO
-
-      ! 2.2 Warm layer; formulation C (Xubin Zeng)
-      !!--------------------------------------------
-      nbi = MAX( 1 , nb_itt_wl)
-
-      IF( nbi > 1 ) THEN
-         !! Itterating for warm-layer solution
-         zdt   = rdt/REAL(nbi)
-         rmult = 1._wp
-      ELSE
-         !! No itteration! Way cheaper !!!
-         !! Assuming balance between the last 2 terms of Eq. 11 (lhs of 11 = 0):
-         !!   => in that case no need for sub-itterations !
-         !! Acceptable in most conditions (UNLESS: sunny + no-wind conditions)
-         !!
-         zdt   = 1._wp
-         rmult = 0._wp
-      END IF
-
-      DO jwl = 1, nbi  ! itteration to solve implicitely equation for warm layer
-
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               ZDSST = pTs(ji,jj) - pSST(ji,jj) - zdT_c(ji,jj)
-
-               !! Buoyancy flux and stability parameter (zdl = -z/L) in water
-               !
-               !! Qt/(rho_w*Cpw):
-               ZSRD = ( pQsw(ji,jj)*ZFI + pQnsol(ji,jj) )/zRhoCp_w
-               !
-               zsgn = 0.5_wp + SIGN(0.5_wp, ZSRD)  ! ZSRD > 0. => 1.  / ZSRD < 0. => 0.
-               ztmp = MAX(ZDSST,0._wp)
-               zdl = (zsgn+1._wp) * ( zus_w2(ji,jj) * SQRT(ztmp/(5._wp*rd0*grav*zalpha_w(ji,jj)/rNu0)) ) & ! (ZDSST > 0.0 .AND. ZSRD < 0.0)
-                  &  +    zsgn    *  ZSRD                                                                  !   otherwize
-               !
-               zus_a = MAX( pustar(ji,jj), 1.E-4_wp )
-               zdL = ZCON3*zalpha_w(ji,jj)*zdL/(zus_a*zus_a*zus_a)
-
-               !! Stability function Phi_t(-z/L) (zdL is -z/L) :
-               zsgn = 0.5_wp + SIGN(0.5_wp, zdL)  ! zdl > 0. => 1.  / zdl < 0. => 0.
-               zdL2 = zdL*zdL
-               ZPHI =    zsgn      * ( 1._wp + (5._wp*zdL + 4._wp*zdL2)/(1._wp + 3._wp*zdL + 0.25_wp*zdL2) ) &  ! (zdL > 0) Takaya et al.
-                  & + (zsgn+1._wp) * ( 1._wp/SQRT(1._wp - 16._wp*(-ABS(zdL))) )        ! (zdl < 0) Eq. 8.136
-               !
-               !! FOR zdL > 0.0, old relations:
-               !         ZPHI = 1.+5._wp*zdL                                ! Eq. 8.136 (Large et al. 1994)
-               !         ZPHI = 1.+5.0*(zdL+zdL**2)/(1.0+3.0*zdL+zdL**2) ! SHEBA, Grachev et al. 2007
-
-               !! Solving 11 by itteration with time step of zdt...
-               ZZ = rmult*1._wp + ZCON4*zdt*zus_w(ji,jj)/ZPHI
-               ZZ = SIGN( MAX(ABS(ZZ) , 1e-4_wp), ZZ )
-               zdT_w(ji,jj) = MAX( 0._wp , (rmult*ZDSST + ZCON5*ZSRD*zdt)/ZZ )
-
-            END DO ! DO ji = 1, jpi
-         END DO ! DO jj = 1, jpj
-
-         ! 3. Apply warm layer and cool skin effects
-         !------------------------------------------
-         pTs(:,:) = pSST(:,:) + zdT_w(:,:) + zdT_c(:,:)
-
-      END DO  ! DO jwl = 1, nbi !: sub-itteration
-
-   END SUBROUTINE CSWL_ECMWF
-
-   !!======================================================================
-END MODULE sbcblk_skin
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare3p6.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare3p6.F90	(revision 11615)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare3p6.F90	(revision 11615)
@@ -0,0 +1,322 @@
+MODULE sbcblk_skin_coare3p6
+   !!======================================================================
+   !!                   ***  MODULE  sbcblk_skin_coare3p6  ***
+   !! Computes:
+   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
+   !!     scheme used at ECMWF
+   !!    Using formulation/param. of COARE 3.6 (Fairall et al., 2019)
+   !!
+   !!   From routine turb_ecmwf maintained and developed in AeroBulk
+   !!            (https://github.com/brodeau/aerobulk)
+   !!
+   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !!----------------------------------------------------------------------
+   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
+   !!               based on the cool-skin/warm-layer scheme used at ECMWF
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+
+   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
+
+   USE lib_mpp         ! distribued memory computing library
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran     ! to use key_nosignedzero
+
+   USE io_ezcdf ! LOLO!
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC :: CS_COARE3P6, WL_COARE3P6
+
+   !! Cool-skin related parameters:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl  !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m]
+
+   !! Warm-layer related parameters:
+   REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp    !: maximum depth of warm layer (adjustable)
+   !
+   REAL(wp), PARAMETER :: rich   = 0.65_wp   !: critical Richardson number
+   !
+   REAL(wp), PARAMETER :: Qabs_thr = 50._wp  !: threshold for heat flux absorbed in WL
+   REAL(wp), PARAMETER :: zfr0   = 0.5_wp    !: initial value of solar flux absorption
+   !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac  ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Qnt_ac    ! time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl     ! depth of warm-layer [m]
+   !!----------------------------------------------------------------------
+CONTAINS
+
+
+   SUBROUTINE CS_COARE3P6( pQsw, pQnsol, pustar, pSST, pQlat,  pdT )
+      !!---------------------------------------------------------------------
+      !!
+      !!  Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019)
+      !!     ------------------------------------------------------------------
+      !!
+      !!  **   INPUT:
+      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
+      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      !!     *pustar*     friction velocity u*                           [m/s]
+      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
+      !!     *pQlat*      surface latent heat flux                       [K]
+      !!
+      !!  **  INPUT/OUTPUT:
+      !!     *pdT*  : as input =>  previous estimate of dT cool-skin
+      !!              as output =>  new estimate of dT cool-skin
+      !!
+      !!------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pSST ! bulk SST [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQlat  ! latent heat flux [W/m^2]
+      !!
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to cool-skin effect
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, &
+         &        zz1, zz2, zus, &
+         &        ztf
+      !!---------------------------------------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
+
+            zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere
+
+            zdelta = delta_vl(ji,jj)   ! using last value of delta
+
+            !! Fraction of the shortwave flux absorbed by the cool-skin sublayer:
+            zfr   = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
+
+            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) )  ! Total cooling at the interface
+
+            ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1
+            !                                 Qt < 0 => warming of the layer => ztf = 0
+
+            !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap):
+            zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap  ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0
+            !! LB: this terms only makes sense if > 0 i.e. in the cooling case
+            !! so similar to what's done in ECMWF:
+            zz1 = MAX(0._wp , zz1)    ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj))
+
+            zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases:
+            zz2 = zus*zus * roadrw
+            zz2 = zz2*zz2
+            zlamb =  6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders)
+
+            ! Updating molecular sublayer thickness (delta):
+            zz2    = rnu0_w/(SQRT(roadrw)*zus)
+            zdelta =      ztf    *          zlamb*zz2   &  ! Eq.12 (when alpha*Qb>0 / cooling of layer)
+               &    + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 )    ! Eq.12 (when alpha*Qb<0 / warming of layer)
+            !LB: changed 0.01 to 0.007
+
+            !! Once again with the new zdelta:
+            zfr   = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq.16 (Fairall al. 1996b)
+            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface
+
+            !! Update!
+            pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment
+            delta_vl(ji,jj) = zdelta
+
+         END DO
+      END DO
+
+   END SUBROUTINE CS_COARE3P6
+
+
+
+
+   SUBROUTINE WL_COARE3P6( kt,  pQsw, pQnsol, pTau, pSST, plon, isd, iwait,  pdT, &
+      &                    Hwl, mask_wl )
+      !!---------------------------------------------------------------------
+      !!
+      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
+      !!     ------------------------------------------------------------------
+      !!
+      !!  **   INPUT:
+      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
+      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      !!     *pTau*       surface wind stress                            [N/m^2]
+      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
+      !!     *plon*       longitude                                      [deg.E]
+      !!     *isd*        current UTC time, counted in second since 00h of the current day
+      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
+      !!
+      !!  **   OUTPUT:
+      !!     *pdT*        dT due to warming at depth of pSST such that SST_actual = pSST + pdT
+      !!---------------------------------------------------------------------
+      !!
+      !!   ** OPTIONAL OUTPUT:
+      !!     *Hwl*        depth of warm layer [m]
+      !!     *mask_wl*    mask for possible existence of a warm-layer (1) or not (0)
+      !!
+      !!------------------------------------------------------------------
+      INTEGER ,                     INTENT(in)  :: kt
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: plon     ! longitude [deg.E]
+      INTEGER ,                     INTENT(in)  :: isd      ! current UTC time, counted in second since 00h of the current day
+      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT      ! dT due to warming at depth of pSST such that pSST_actual = pSST + pdT
+      !!
+      REAL(wp),   DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl     ! depth of warm layer [m]
+      INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0)
+      !
+      !
+      INTEGER :: ji,jj
+      !
+      REAL(wp) :: zdT_wl, zQabs, zfr, zdz
+      REAL(wp) :: zqac, ztac
+      REAL(wp) :: zalpha_w, zcd1, zcd2, flg
+      CHARACTER(len=128) :: cf_tmp
+      !!---------------------------------------------------------------------
+
+      REAL(wp) :: rlag_gw_h, &  ! local solar time lag in hours   / Greenwich meridian (lon==0) => ex: ~ -10.47 hours for Hawai
+         &        rhr_sol       ! local solar time in hours since midnight
+
+      INTEGER  :: ilag_gw_s, &  ! local solar time LAG in seconds / Greenwich meridian (lon==0) => ex: ~ INT( -10.47*3600. ) seconds for Hawai
+         &        isd_sol,   &  ! local solar time in number of seconds since local solar midnight
+         &        jl
+
+      !! INITIALIZATION:
+      pdT(:,:) = 0._wp    ! dT initially set to 0._wp
+      zdT_wl = 0._wp       ! total warming (amplitude) in warm layer
+      zQabs  = 0._wp       ! total heat absorped in warm layer
+      zfr    = zfr0        ! initial value of solar flux absorption
+      ztac   = 0._wp
+      zqac   = 0._wp
+      IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            zdz   = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer
+
+            !! Need to know the local solar time from longitude and isd:
+            rlag_gw_h = -1._wp * MODULO( ( 360._wp - MODULO(plon(ji,jj),360._wp) ) / 15._wp , 24._wp )
+            rlag_gw_h = -1._wp * SIGN( MIN(ABS(rlag_gw_h) , ABS(MODULO(rlag_gw_h,24._wp))), rlag_gw_h + 12._wp )
+            ilag_gw_s = INT( rlag_gw_h*3600._wp )
+            isd_sol = MODULO( isd + ilag_gw_s , 24*3600 )
+            rhr_sol = REAL( isd_sol , wp) / 3600._wp
+            !LOLO!:
+            !pdT(ji,jj) = REAL( isd_sol , wp) / 3600._wp !LOLO!
+            !END DO
+            !END DO
+            !CALL DUMP_FIELD(REAL( pdT , 4), 'isd_sol.nc', 'solar_hour' )
+            !STOP'LOLO!'
+            !DO jj = 1, jpj
+            !   DO ji = 1, jpi
+            !LOLO.
+
+            !PRINT *, ' Lag in hours / Greenwich for local solar time =', rlag_gw_h
+            !PRINT *, '     UTC     time in seconds:', isd
+            !PRINT *, ' Local solar time in seconds:', isd_sol
+            !PRINT *, '     UTC     time in hours:',   REAL(isd    ,wp)/3600._wp
+            !PRINT *, '  [WL_COARE3P6] Local solar time in hours:',   REAL(isd_sol,wp)/3600._wp
+
+            !*****  variables for warm layer  ***
+            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
+
+            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w))        !mess-o-constants 1
+            zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
+
+            !********************************************************
+            !****  Compute apply warm layer  correction *************
+            !********************************************************
+            
+            !IF (isd_sol <= rdt ) THEN    !re-zero at midnight ! LOLO improve: risky if real midnight (00:00:00) is not a time in vtime...
+            IF ( (rhr_sol > 23.5_wp).OR.(rhr_sol < 4._wp) ) THEN
+               !PRINT *, '  [WL_COARE3P6] MIDNIGHT RESET !!!!, isd_sol =>', isd_sol
+               zdz           = H_wl_max
+               Tau_ac(ji,jj) = 0._wp
+               Qnt_ac(ji,jj) = 0._wp
+            END IF
+
+            IF ( rhr_sol > 5._wp ) THEN  ! ( 5am)
+
+               !PRINT *, '  [WL_COARE3P6] WE DO WL !'
+               !PRINT *, '  [WL_COARE3P6] isd_sol, pTau, pSST, pdT =', isd_sol, REAL(pTau(ji,jj),4), REAL(pSST(ji,jj),4), REAL(pdT(ji,jj),4)
+
+               !************************************
+               !****   set warm layer constants  ***
+               !************************************
+
+               zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
+
+               !PRINT *, '  [WL_COARE3P6] rdt,  pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4)
+
+               IF ( zQabs >= Qabs_thr ) THEN         ! Check for threshold
+
+                  !PRINT *, '  [WL_COARE3P6] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4)
+
+                  !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! momentum integral
+                  ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
+
+                  IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN         !check threshold for warm layer existence
+                     !******************************************
+                     ! Compute the absorption profile
+                     !******************************************
+                     DO jl = 1, 5                           !loop 5 times for zfr
+                        zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) &
+                           &        + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz
+                        zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
+                        IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
+                     END DO
+
+                  ELSE
+                     !***********************
+                     ! Warm layer wiped out
+                     !***********************
+                     zfr  = 0.75
+                     zdz  = H_wl_max
+                     zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
+
+                  END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp )
+
+                  IF ( zqac > 1._wp ) zdT_wl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp)  !! => IF(zqac>0._wp): zdT_wl=zcd2*zqac**1.5/ztac ; ELSE: zdT_wl=0. / ! normally: zqac > 0 !
+                  
+                  ! Warm layer correction
+                  flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = zdT_wl) | 0 when gdept_1d(1)<zdz (pdT(ji,jj) = zdT_wl*gdept_1d(1)/zdz)
+                  pdT(ji,jj) = zdT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
+                  
+               END IF ! IF ( zQabs >= Qabs_thr )
+
+            END IF ! IF ( isd_sol >= 21600 ) THEN  ! (21600 == 6am)
+
+            IF ( iwait == 0 ) THEN
+               IF ( (zQabs >= Qabs_thr).AND.(rhr_sol >= 5._wp) ) THEN
+                  !PRINT *, '  [WL_COARE3P6] WE UPDATE ACCUMULATED FLUXES !!!'
+                  Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
+                  Tau_ac(ji,jj) = ztac !
+                  IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1
+               END IF
+            END IF
+
+            H_wl(ji,jj) = zdz
+
+            IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj)
+
+         END DO
+      END DO
+
+      !WRITE(cf_tmp,'("Qnt_ac_",i5.5,".nc")') kt
+      !CALL DUMP_FIELD(REAL( Qnt_ac , 4), TRIM(cf_tmp), 'Qnt_ac')
+      
+      
+   END SUBROUTINE WL_COARE3P6
+
+
+   !!======================================================================
+END MODULE sbcblk_skin_coare3p6
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11615)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11615)
@@ -0,0 +1,223 @@
+MODULE sbcblk_skin_ecmwf
+   !!======================================================================
+   !!                   ***  MODULE  sbcblk_skin_ecmwf  ***
+   !! Computes:
+   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
+   !!     scheme used at ECMWF
+   !!    Using formulation/param. of IFS of ECMWF (cycle 40r1)
+   !!   based on IFS doc (avaible online on the ECMWF's website)
+   !!
+   !!   From routine turb_ecmwf maintained and developed in AeroBulk
+   !!            (https://github.com/brodeau/aerobulk)
+   !!
+   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !!----------------------------------------------------------------------
+   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
+   !!               based on the cool-skin/warm-layer scheme used at ECMWF
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+
+   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
+
+   USE lib_mpp         ! distribued memory computing library
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran     ! to use key_nosignedzero
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC :: CS_ECMWF, WL_ECMWF
+
+   !! Cool-skin related parameters:
+   ! ...
+
+   !! Warm-layer related parameters:
+   REAL(wp), PARAMETER :: rd0  = 3.        !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005)
+   REAL(wp), PARAMETER :: zRhoCp_w = rho0_w*rCp0_w
+   REAL(wp), PARAMETER :: rNu0 = 1.0       !:  be closer to COARE3p6 ???!LOLO
+   !REAL(wp), PARAMETER :: rNu0 = 0.5       !: Nu (exponent of temperature profile) Eq.11
+   !                                       !: (Zeng & Beljaars 2005) !: set to 0.5 instead of
+   !                                       !: 0.3 to respect a warming of +3 K in calm
+   !                                       !: condition for the insolation peak of +1000W/m^2
+   !!----------------------------------------------------------------------
+CONTAINS
+
+
+   SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST,  pdT )
+      !!---------------------------------------------------------------------
+      !!
+      !!  Cool-Skin scheme according to Fairall et al. 1996
+      !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
+      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
+      !!     as parameterized in IFS Cy45r1 / E.C.M.W.F.
+      !!     ------------------------------------------------------------------
+      !!
+      !!  **   INPUT:
+      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
+      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      !!     *pustar*     friction velocity u*                           [m/s]
+      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
+      !!
+      !!  **  INPUT/OUTPUT:
+      !!     *pdT*  : as input =>  previous estimate of dT cool-skin
+      !!              as output =>  new estimate of dT cool-skin
+      !!
+      !!------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pSST ! bulk SST [K]
+      !!
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to cool-skin effect
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, &
+         & zusw, zusw2
+      !!---------------------------------------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
+
+            zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere
+
+            zusw  = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)    ! u* in the water
+            zusw2 = zusw*zusw
+
+            zlamb = 6._wp*( 1._wp + (zQnsol*zalpha_w*rcst_cs/(zusw2*zusw2 ))**0.75 )**(-1./3.)   ! w.r.t COARE 3p6 => seems to ommit absorbed zfr*Qsw (Qnet i.o. Qnsol) and effect of evap...
+            !                                                                                  ! so zlamb not implicit in terms of zdelta (zfr(delta)), so no need to have last guess of delta as in COARE 3.6...
+            zdelta = zlamb*rnu0_w/zusw
+
+            zfr   = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq. 8.131 / IFS cy40r1, doc, Part IV,
+            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface
+
+            !! Update!
+            pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment
+
+         END DO
+      END DO
+
+   END SUBROUTINE CS_ECMWF
+
+
+
+   SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pdT )
+      !!---------------------------------------------------------------------
+      !!
+      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
+      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
+      !!
+      !!    As included in IFS Cy45r1   /  E.C.M.W.F.
+      !!     ------------------------------------------------------------------
+      !!
+      !!  **   INPUT:
+      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
+      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      !!     *pustar*     friction velocity u*                           [m/s]
+      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
+      !!
+      !!   **  INPUT/OUTPUT:
+      !!     *pdT*  : as input =>  previous estimate of dT warm-layer
+      !!             as output =>  new estimate of dT warm-layer
+      !!
+      !!------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar   ! friction velocity [m/s]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
+      !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! difference between "almost surface (right below viscous layer) and depth of bulk SST
+      !
+      INTEGER :: ji,jj
+      !
+      REAL(wp) :: &
+         & zdz,    & !: thickness of the warm-layer [m]
+         & zalpha_w, & !: thermal expansion coefficient of sea-water
+         & ZSRD,   &
+         & dT_wl,   & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL
+         & zfr,zdL,zdL2, ztmp, &
+         & ZPHI, &
+         & zus_a, zusw, zusw2, &
+         & flg, zQabs, ZL1, ZL2
+      !!---------------------------------------------------------------------
+
+      IF (lwp) WRITE(numout,*) 'LOLO: WL_ECMWF@sbcblk_skin_ecmwf.F90: depth for SST bulk takent at', gdept_1d(1), ' m'
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            zdz = rd0 ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !)
+
+            ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz)
+            ! pdT         "                          "                                    and depth of bulk SST (here gdept_1d(1))!
+            !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced!
+            !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof dT_wl !
+            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$
+            dT_wl = pdT(ji,jj) / ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
+            !PRINT *, 'LOLO/mod_wl_ecmwf.f90: dT_wl2=', dT_wl
+            !PRINT *, ''
+
+            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
+
+
+            ! *** zfr = Fraction of solar radiation absorbed in warm layer (-)
+            zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*zdz) - 0.27_wp*EXP(-2.8_wp*zdz) - 0.45_wp*EXP(-0.07_wp*zdz)  !: Eq. 8.157
+
+            zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
+
+            zusw  = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)    ! u* in the water
+            zusw2 = zusw*zusw
+
+
+            !! *** 1st rhs term in eq. 8.156 (IFS doc Cy45r1):
+            ZL1 = zQabs / ( zdz * zRhoCp_w * rNu0 ) * (rNu0 + 1._wp)
+
+
+            !! Buoyancy flux and stability parameter (zdl = -z/L) in water
+            ZSRD = zQabs/zRhoCp_w
+            !
+            flg = 0.5_wp + SIGN(0.5_wp, ZSRD)  ! ZSRD > 0. => 1.  / ZSRD < 0. => 0.
+            ztmp = MAX(dT_wl,0._wp)
+            zdl = (flg+1._wp) * ( zusw2 * SQRT(ztmp/(5._wp*zdz*grav*zalpha_w/rNu0)) ) & ! (dT_wl > 0.0 .AND. ZSRD < 0.0)
+               &  +    flg    *  ZSRD                                                                  !   otherwize
+            !
+            zus_a = MAX( pustar(ji,jj), 1.E-4_wp )
+            zdL = zdz*vkarmn*grav/(roadrw)**1.5_wp*zalpha_w*zdL/(zus_a*zus_a*zus_a)
+
+            !! Stability function Phi_t(-z/L) (zdL is -z/L) :
+            flg = 0.5_wp + SIGN(0.5_wp, zdL)  ! zdl > 0. => 1.  / zdl < 0. => 0.
+            zdL2 = zdL*zdL
+            ZPHI =    flg      * ( 1._wp + (5._wp*zdL + 4._wp*zdL2)/(1._wp + 3._wp*zdL + 0.25_wp*zdL2) ) &  ! (zdL > 0) Takaya et al.
+               & + (flg+1._wp) * ( 1._wp/SQRT(1._wp - 16._wp*(-ABS(zdL))) )        ! (zdl < 0) Eq. 8.136
+            !! FOR zdL > 0.0, old relations:
+            !         ZPHI = 1.+5._wp*zdL                                ! Eq. 8.136 (Large et al. 1994)
+            !         ZPHI = 1.+5.0*(zdL+zdL**2)/(1.0+3.0*zdL+zdL**2) ! SHEBA, Grachev et al. 2007
+
+            !! *** 2nd rhs term in eq. 8.156 (IFS doc Cy45r1):
+            ZL2 = - (rNu0 + 1._wp) * vkarmn * zusw / ( zdz * ZPHI )
+
+            ! Forward time / explicit solving of eq. 8.156 (IFS doc Cy45r1): (f_n+1 == pdT(ji,jj) ; f_n == dT_wl)
+            dT_wl = MAX ( dT_wl + rdt*ZL1 + rdt*ZL2*dT_wl , 0._wp )
+
+            ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz)
+            !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced!
+
+            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$
+            pdT(ji,jj) = dT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
+
+         END DO
+      END DO
+
+   END SUBROUTINE WL_ECMWF
+
+   !!======================================================================
+END MODULE sbcblk_skin_ecmwf
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/demo_cfgs.txt
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/demo_cfgs.txt	(revision 11614)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/demo_cfgs.txt	(revision 11615)
@@ -9,2 +9,3 @@
 WAD OCE
 BENCH OCE ICE TOP
+STATION_TASF OCE
