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 11265)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90	(revision 11266)
@@ -153,5 +153,5 @@
    !!                     Cool-skin/Warm-layer
    !!----------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature (used if ln_skin==.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 11265)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11266)
@@ -346,6 +346,7 @@
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
       !
-      !                                            ! compute the surface ocean fluxes using bulk formulea
-      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
+      IF( kt == nit000 ) tsk(:,:) = sst_m(:,:)*tmask(:,:,1)  ! no previous estimate of skin temperature => using bulk SST
+      !
+      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) !  compute the surface ocean fluxes using bulk formulea
 
 #if defined key_cice
@@ -493,11 +494,13 @@
             CALL turb_coare   ( 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))
+               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
+               &                Tsk_b=tsk(:,:) )
 
          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))
-
+               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
+               &                Tsk_b=tsk(:,:) )
+            
          CASE DEFAULT
             CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' )
@@ -515,5 +518,5 @@
          END WHERE
 
-         !LB: Update of tsk, the official array for skin temperature
+         !LB: Update of tsk, the "official" array for skin temperature
          tsk(:,:) = zst(:,:)
 
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11265)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11266)
@@ -57,5 +57,5 @@
       &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
       &                   Cdn, Chn, Cen,                      &
-      &                   Qsw, rad_lw, slp                   )
+      &                   Qsw, rad_lw, slp, Tsk_b             )
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_coare  ***
@@ -81,5 +81,8 @@
       !! INPUT/OUTPUT:
       !! -------------
-      !!    *  T_s  : SST or skin temperature                                 [K]
+      !!    *  T_s  : always "bulk SST" as input                              [K]
+      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
+      !!              -> skin temperature as output if CSWL used              [K]
+      !!
       !!    *  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)
@@ -91,4 +94,5 @@
       !!    *  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]
       !!
       !! OUTPUT :
@@ -122,4 +126,5 @@
       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]
       !
       INTEGER :: j_itt
@@ -152,7 +157,12 @@
 
       !! Initialization for cool skin:
+      zsst   = T_s    ! save the bulk SST
       IF( l_use_skin ) THEN
-         zsst   = T_s    ! save the bulk SST
-         T_s    = T_s - 0.25                      ! First guess of correction
+         ! 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
@@ -202,5 +212,5 @@
       q_star = dq_zu*ztmp0
 
-      ! What's need to be done if zt /= zu:
+      ! What needs to be done if zt /= zu:
       IF( .NOT. l_zt_equal_zu ) THEN
          !! First update of values at zu (or zt for wind)
@@ -275,6 +285,7 @@
             ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
             ztmp2 = T_s*T_s
-            ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux
-               &       +    (rad_lw - emiss_w*stefan*ztmp2*ztmp2)                  ! Net longwave flux
+            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
+               &       +      rad_lw - emiss_w*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!)
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 11265)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11266)
@@ -63,5 +63,5 @@
       &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
       &                   Cdn, Chn, Cen,                      &
-      &                   Qsw, rad_lw, slp                   )
+      &                   Qsw, rad_lw, slp, Tsk_b             )
       !!----------------------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_ecmwf  ***
@@ -87,5 +87,8 @@
       !! INPUT/OUTPUT:
       !! -------------
-      !!    *  T_s  : SST or skin temperature                                 [K]
+      !!    *  T_s  : always "bulk SST" as input                              [K]
+      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
+      !!              -> skin temperature as output if CSWL used              [K]
+      !!
       !!    *  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)
@@ -97,4 +100,5 @@
       !!    *  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]
       !!
       !! OUTPUT :
@@ -128,4 +132,5 @@
       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]
       !
       INTEGER :: j_itt
@@ -156,10 +161,15 @@
       !
       l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01 )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+      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
-         zsst   = T_s    ! save the bulk SST
-         T_s    = T_s - 0.25                      ! First guess of correction
+         ! 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
@@ -209,5 +219,5 @@
       q_star = dq_zu*ztmp0
 
-      ! What's need to be done if zt /= zu:
+      ! What needs to be done if zt /= zu:
       IF( .NOT. l_zt_equal_zu ) THEN
          !! First update of values at zu (or zt for wind)
@@ -304,6 +314,7 @@
             ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
             ztmp2 = T_s*T_s
-            ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux
-               &       +    (rad_lw - emiss_w*stefan*ztmp2*ztmp2)                  ! Net longwave flux
+            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
+               &       +      rad_lw - emiss_w*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 )
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 11265)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11266)
@@ -219,5 +219,10 @@
             zqa = (1._wp + rctv0*pqa(ji,jj))
             !
-            One_on_L(ji,jj) = grav*vkarmn*(pts(ji,jj) + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
+            ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
+            !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
+            !                      or
+            !  b/  -u* [ theta*              + 0.61 theta q* ]
+            !
+            One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
                &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
             !
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 11265)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin.F90	(revision 11266)
@@ -34,9 +34,9 @@
    IMPLICIT NONE
    PRIVATE
-   
-   PUBLIC :: CSWL_ECMWF   ! called by sbcblk_algo_*.F90
+
+   PUBLIC :: CSWL_ECMWF   ! called by sbcblk_algo_coare.F90 and sbcblk_algo_ecmwf.F90
 
    !! Cool-Skin / Warm-Layer related parameters:
-   REAL(wp), PARAMETER :: rdt0 = 3600.*1.5 !: time step
+   !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
@@ -162,5 +162,5 @@
       IF( nbi > 1 ) THEN
          !! Itterating for warm-layer solution
-         zdt   = rdt0/REAL(nbi)
+         zdt   = rdt/REAL(nbi)
          rmult = 1._wp
       ELSE
