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 11110)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/DOM/phycst.F90	(revision 11111)
@@ -50,5 +50,5 @@
    REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice (for pisces)          [psu]
    REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea (for pisces and isf)  [psu]
-   REAL(wp), PUBLIC ::   rLevap   =    2.5e+6_wp     !: latent heat of evaporation (water)
+   REAL(wp), PUBLIC ::   rLevap   =    2.46e+6_wp    !: latent heat of vaporization for sea-water   [J/kg] #LB
    REAL(wp), PUBLIC ::   vkarmn   =    0.4_wp        !: von Karman constant
    REAL(wp), PUBLIC ::   stefan   =    5.67e-8_wp    !: Stefan-Boltzmann constant 
@@ -66,4 +66,33 @@
    REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos
    REAL(wp), PUBLIC ::   r1_rcpi                     !: 1 / rcpi
+
+
+   !#LB:
+   !! Mainly used in OCE/SBC (removed from sbcblk.F90)
+   REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp   !: Specic heat of dry air, constant pressure      [J/K/kg]
+   REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp   !: Specic heat of water vapor, constant pressure  [J/K/kg]
+   REAL(wp), PARAMETER, PUBLIC :: R_dry   = 287.05_wp   !: Specific gas constant for dry air              [J/K/kg]
+   REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg]
+   REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
+   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+   REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...)
+   REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp   !: transfer coefficient over ice
+   REAL(wp), PARAMETER, PUBLIC :: albo    = 0.066_wp    !: ocean albedo assumed to be constant
+   !
+   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 :: 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]
+   REAL(wp), PARAMETER, PUBLIC :: rk0_w   = 0.6_wp      !: thermal conductivity of water (at 20C)          [W/m/K]
+   !
+   REAL(wp), PARAMETER, PUBLIC :: emiss_w = 1._wp       !: Surface emissivity (black-body long-wave radiation) of sea-water []
+   !                                                    !: Theoretically close to 0.97! Yet, taken equal as 1 to account for
+   !                                                    !: 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  
+   !#LB.
+
+
+   
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
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 11110)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90	(revision 11111)
@@ -149,4 +149,10 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-]
 
+   !!----------------------------------------------------------------------
+   !!                     Cool-skin/Warm-layer
+   !!----------------------------------------------------------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius]
+
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
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 11110)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11111)
@@ -15,5 +15,5 @@
    !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
    !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
-   !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
+   !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
    !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
    !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle) 
@@ -24,8 +24,4 @@
    !!   sbc_blk       : bulk formulation as ocean surface boundary condition
    !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
-   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
-   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
-   !!   q_sat         : saturation humidity as a function of SLP and temperature
-   !!   L_vap         : latent heat of vaporization of water as a function of temperature
    !!             sea-ice case only : 
    !!   blk_ice_tau   : provide the air-ice stress
@@ -60,4 +56,7 @@
    USE prtctl         ! Print control
 
+   USE sbcblk_phymbl  !LB: all thermodynamics functions, rho_air, q_sat, etc...
+
+   
    IMPLICIT NONE
    PRIVATE
@@ -70,13 +69,4 @@
    PUBLIC   blk_ice_qcn   ! routine called in icesbc
 #endif 
-
-!!Lolo: should ultimately be moved in the module with all physical constants ?
-!!gm  : In principle, yes.
-   REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg]
-   REAL(wp), PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg]
-   REAL(wp), PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
-   REAL(wp), PARAMETER ::   R_vap = 461.495_wp    !: Specific gas constant for water vapor          [J/K/kg]
-   REAL(wp), PARAMETER ::   reps0 = R_dry/R_vap   !: ratio of gas constant for dry air and water vapor => ~ 0.622
-   REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
 
    INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
@@ -94,11 +84,4 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
 
-   !                                             !!! Bulk parameters
-   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...)
-   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation
-   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant
-   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice
-   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant
-   !
    !                           !!* Namelist namsbc_blk : bulk parameters
    LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
@@ -124,4 +107,9 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
 
+   !LB:
+   LOGICAL  ::   ln_skin        ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
+   !LB: => sbc_oce.F90: REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk                       ! sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius]
+   !LB.
+   
    INTEGER  ::   nblk           ! choice of the bulk algorithm
    !                            ! associated indices:
@@ -150,4 +138,17 @@
       IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
    END FUNCTION sbc_blk_alloc
+
+   !LB:
+   INTEGER FUNCTION sbc_blk_cswl_alloc()
+      !!-------------------------------------------------------------------
+      !!             ***  ROUTINE sbc_blk_cswl_alloc ***
+      !!-------------------------------------------------------------------
+      !PRINT *, '*** LB: allocating tsk!'
+      ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc )
+      !PRINT *, '*** LB: done!'
+      CALL mpp_sum ( 'sbcblk', sbc_blk_cswl_alloc )
+      IF( sbc_blk_cswl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_cswl_alloc: failed to allocate arrays' )
+   END FUNCTION sbc_blk_cswl_alloc
+   !LB.
 
 
@@ -173,5 +174,6 @@
          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
-         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
+         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
+         &                 ln_skin                                                        ! cool-skin / warm-layer param. !LB
       !!---------------------------------------------------------------------
       !
@@ -198,4 +200,13 @@
       IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF
       !
+      !LB:
+      !                             !** initialization of the cool-skin / warm-layer parametrization
+      IF( ln_skin ) 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.
+         IF( sbc_blk_cswl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
+      END IF
+      !LB.
+         
       IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
       !
@@ -279,4 +290,7 @@
          END SELECT
          !
+         WRITE(numout,*)
+         WRITE(numout,*) '      use cool-skin/warm-layer parameterization (SSST)     ln_skin     = ', ln_skin !LB
+         !
       ENDIF
       !
@@ -413,5 +427,5 @@
 
       ! ----------------------------------------------------------------------------- !
-      !      I   Radiative FLUXES                                                     !
+      !      I   Solar FLUX                                                           !
       ! ----------------------------------------------------------------------------- !
 
@@ -422,5 +436,4 @@
       ENDIF
 
-      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
 
       ! ----------------------------------------------------------------------------- !
@@ -429,5 +442,5 @@
 
       ! ... specific humidity at SST and IST tmask(
-      zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
+      zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
       !!
       !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
@@ -444,9 +457,32 @@
       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+
+      !LB: Skin!!!
+      CASE( np_ECMWF     )
+         IF ( ln_skin ) THEN            
+            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL optional arrays!!!'
+            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => BEFORE ZST(40:50,30) =', zst(40:50,30)
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF         
+               &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
+               &                                           Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1))
+            !LB: "zst" and "zsq" have been updated with skin temp. !!! (from bulk SST)...
+            zst(:,:) = zst(:,:)*tmask(:,:,1)
+            zsq(:,:) = zsq(:,:)*tmask(:,:,1)
+            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => AFTER ZST(40:50,30) =', zst(40:50,30)
+         ELSE
+            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF         
+               &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         END IF
+         !LB.
+         
       CASE DEFAULT
          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
       END SELECT
+
+
+      !LB: cool-skin/warm-layer:
+      IF ( ln_skin ) tsk(:,:) = zst(:,:)
+
 
       !                          ! Compute true air density :
@@ -507,4 +543,12 @@
 
 
+      ! ----------------------------------------------------------------------------- !
+      !     III    Net longwave radiative FLUX                                        !
+      ! ----------------------------------------------------------------------------- !
+
+      !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin==.TRUE.)
+      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - emiss_w * stefan * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
+
+
       IF(ln_ctl) THEN
          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
@@ -519,5 +563,5 @@
 
       ! ----------------------------------------------------------------------------- !
-      !     III    Total FLUXES                                                       !
+      !     IV    Total FLUXES                                                       !
       ! ----------------------------------------------------------------------------- !
       !
@@ -527,5 +571,5 @@
       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
-         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
+         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !?
          &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
          &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
@@ -539,5 +583,6 @@
 #endif
       !
-      IF ( nn_ice == 0 ) THEN
+      !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN
+         CALL iom_put( "rho_air" ,   zrhoa )                ! output air density (kg/m^3) !#LB
          CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
          CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
@@ -551,5 +596,9 @@
          CALL iom_put( 'snowpre', sprecip )                 ! Snow
          CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
-      ENDIF
+         IF ( ln_skin ) 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...
+         END IF
+      !!#LB. ENDIF
       !
       IF(ln_ctl) THEN
@@ -564,121 +613,4 @@
 
 
-
-   FUNCTION rho_air( ptak, pqa, pslp )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION rho_air  ***
-      !!
-      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
-      !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
-      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
-
-
-   FUNCTION cp_air( pqa )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION cp_air  ***
-      !!
-      !! ** Purpose : Compute specific heat (Cp) of moist air
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!-------------------------------------------------------------------------------
-      REAL(wp), 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 = Cp_dry + Cp_vap * pqa
-      !
-   END FUNCTION cp_air
-
-
-   FUNCTION q_sat( ptak, pslp )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION q_sat  ***
-      !!
-      !! ** Purpose : Specific humidity at saturation in [kg/kg]
-      !!              Based on accurate estimate of "e_sat"
-      !!              aka saturation water vapor (Goff, 1957)
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::   ze_sat, ztmp   ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            ztmp = rt0 / ptak(ji,jj)
-            !
-            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
-            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
-               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
-               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
-               !
-            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION q_sat
-
-
-   FUNCTION gamma_moist( ptak, pqa )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION gamma_moist  ***
-      !!
-      !! ** Purpose : Compute the moist adiabatic lapse-rate.
-      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
-      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) :: zrv, ziRT        ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
-            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
-            gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )
-         END DO
-      END DO
-      !
-   END FUNCTION gamma_moist
-
-
-   FUNCTION L_vap( psst )
-      !!---------------------------------------------------------------------------------
-      !!                           ***  FUNCTION L_vap  ***
-      !!
-      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! 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
 
 #if defined key_si3
@@ -710,7 +642,7 @@
       !
       ! set transfer coefficients to default sea-ice values
-      Cd_atm(:,:) = Cd_ice
-      Ch_atm(:,:) = Cd_ice
-      Ce_atm(:,:) = Cd_ice
+      Cd_atm(:,:) = rCd_ice
+      Ch_atm(:,:) = rCd_ice
+      Ce_atm(:,:) = rCd_ice
 
       wndm_ice(:,:) = 0._wp      !!gm brutal....
@@ -737,5 +669,5 @@
       ENDIF
 
-!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
+!!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
 !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
 
@@ -792,5 +724,5 @@
       REAL(wp) ::   zst3                     ! local variable
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
-      REAL(wp) ::   zztmp, z1_rLsub           !   -      -
+      REAL(wp) ::   zztmp, z1_rLsub          !   -      -
       REAL(wp) ::   zfr1, zfr2               ! local variables
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
@@ -803,6 +735,6 @@
       !!---------------------------------------------------------------------
       !
-      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars
-      zcoef_dqla = -Ls * 11637800. * (-5897.8)
+      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars
+      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
       !
       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
@@ -824,5 +756,5 @@
                qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
                ! Long  Wave (lw)
-               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
+               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ! lw sensitivity
                z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
@@ -834,7 +766,7 @@
                ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
                ! Sensible Heat
-               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
+               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
                ! Latent Heat
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
                   &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
                ! Latent heat sensitivity for ice (Dqla/Dt)
@@ -847,5 +779,5 @@
 
                ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
-               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
+               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
 
                ! ----------------------------!
@@ -1064,9 +996,9 @@
       ! generic drag over a cell partly covered by ice
       !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
-      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
+      !!   &      rCd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
       !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
 
       ! ice-atm drag
-      Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag
+      Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag
          &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
       
@@ -1141,6 +1073,6 @@
       ! Atmospheric and Surface Variables
       zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
-      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
-      zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
+      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
+      zqi_sat(:,:) =                  q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt)
       !
       DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
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 11110)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11111)
@@ -1,8 +1,6 @@
 MODULE sbcblk_algo_coare
    !!======================================================================
-   !!                       ***  MODULE  sbcblk_algo_coare  ***
-   !! Computes turbulent components of surface fluxes
-   !!         according to Fairall et al. 2003 (COARE v3)
-   !!
+   !!                   ***  MODULE  sbcblk_algo_coare  ***
+   !! Computes:
    !!   * bulk transfer coefficients C_D, C_E and C_H
    !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
@@ -10,14 +8,12 @@
    !!   => all these are used in bulk formulas in sbcblk.F90
    !!
-   !!    Using the bulk formulation/param. of COARE v3, Fairall et al. 2003
-   !!
+   !!    Using the bulk formulation/param. of COARE v3, Fairall et al., 2003
    !!
    !!       Routine turb_coare maintained and developed in AeroBulk
-   !!                     (http://aerobulk.sourceforge.net/)
+   !!                     (https://github.com/brodeau/aerobulk)
    !!
-   !!            Author: Laurent Brodeau, 2016, brodeau@gmail.com
-   !!
-   !!======================================================================
-   !! History :  3.6  !  2016-02  (L.Brodeau)   Original code
+   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !!----------------------------------------------------------------------
+   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
    !!----------------------------------------------------------------------
 
@@ -27,13 +23,13 @@
    !!                   returns the effective bulk wind speed at 10m
    !!----------------------------------------------------------------------
-   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 sbcwave, ONLY  :  cdn_wave ! wave module
+   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 sbcwave, ONLY   :  cdn_wave ! wave module
 #if defined key_si3 || defined key_cice
-   USE sbc_ice        ! Surface boundary condition: ice fields
+   USE sbc_ice         ! Surface boundary condition: ice fields
 #endif
-   !
+   USE sbcblk_phymbl  !LB: all thermodynamics functions, rho_air, q_sat, etc...
    USE in_out_manager ! I/O manager
    USE iom            ! I/O manager library
@@ -50,5 +46,6 @@
    REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
    REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
-   REAL(wp), PARAMETER ::   rctv0   =   0.608_wp   ! constant to obtain virtual temperature...
+
+   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
 
    !!----------------------------------------------------------------------
@@ -61,12 +58,9 @@
       !!                      ***  ROUTINE  turb_coare  ***
       !!
-      !!            2015: L. Brodeau (brodeau@gmail.com)
-      !!
       !! ** Purpose :   Computes turbulent transfert coefficients of surface
       !!                fluxes according to Fairall et al. (2003)
       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
-      !!
-      !! ** Method : Monin Obukhov Similarity Theory
-      !!----------------------------------------------------------------------
+      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
+      !!
       !!
       !! INPUT :
@@ -88,6 +82,9 @@
       !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
       !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
-      !!    *  U_blk  : bulk wind at 10m                                      [m/s]
-      !!----------------------------------------------------------------------
+      !!    *  U_blk  : bulk wind speed at 10m                                [m/s]
+      !!
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
+      !!----------------------------------------------------------------------------------
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
       REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
@@ -106,6 +103,5 @@
       !
       INTEGER :: j_itt
-      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
-      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
+      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
 
       REAL(wp), DIMENSION(jpi,jpj) ::  &
@@ -125,12 +121,12 @@
 
       !! First guess of temperature and humidity at height zu:
-      t_zu = MAX(t_zt , 0.0)    ! who knows what's given on masked-continental regions...
-      q_zu = MAX(q_zt , 1.e-6)  !               "
+      t_zu = MAX( t_zt , 199.0_wp )   ! who knows what's given on masked-continental regions...
+      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
 
       !! Pot. temp. difference (and we don't want it to be 0!)
-      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
-      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
-
-      znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K)
+      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
+      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
+
+      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
 
       ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
@@ -144,21 +140,24 @@
 
       z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star
-      z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t !
+      z0     = MIN(ABS(z0), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
+      z0t    = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
+      z0t    = MIN(ABS(z0t), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
       ztmp2  = vkarmn/ztmp0
       Cd     = ztmp2*ztmp2    ! first guess of Cd
 
-      ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd
-
-      !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
-      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number
-      ztmp1 = 0.5 + sign(0.5 , ztmp2)
+      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
+
+      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
+
+      !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2):
+      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
       ztmp0 = ztmp0*ztmp2
-      !!             Ribu < 0                                 Ribu > 0   Beta = 1.25
-      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &
-         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))
+      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  Ribu < 0
+         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  Ribu > 0
+      !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
 
       !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
-      ztmp0  =  vkarmn/(LOG(zu*z0t) - psi_h_coare(zeta_u))
+      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
 
       u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
@@ -168,17 +167,14 @@
       ! What's need to be done if zt /= zu:
       IF( .NOT. l_zt_equal_zu ) THEN
-
+         !! First update of values at zu (or zt for wind)
          zeta_t = zt*zeta_u/zu
-
-         !! First update of values at zu (or zt for wind)
          ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
-         ztmp1 = log(zt/zu) + ztmp0
+         ztmp1 = LOG(zt/zu) + ztmp0
          t_zu = t_zt - t_star/vkarmn*ztmp1
          q_zu = q_zt - q_star/vkarmn*ztmp1
-         q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity :
-
-         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
-         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
-
+         q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
+         !
+         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
+         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
       END IF
 
@@ -188,4 +184,5 @@
          !!Inverse of Monin-Obukov length (1/L) :
          ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
+         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO
 
          ztmp1 = u_star*u_star   ! u*^2
@@ -193,7 +190,7 @@
          !! Update wind at 10m 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.))**(2./3.)   ! => ztmp2 == Ug^2
+         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.)   ! => 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)        ! include gustiness in bulk wind speed
+         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.
 
@@ -203,20 +200,22 @@
          !! Roughness lengthes z0, z0t (z0q = z0t) :
          z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
-         ztmp1 = z0*u_star/znu_a                             ! Re_r: roughness Reynolds number
-         z0t  = min( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )         ! Scalar roughness for both theta and q (eq.28)
+         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)
 
          !! Stability parameters:
-         zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u )
+         zeta_u = zu*ztmp0
+         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
          IF( .NOT. l_zt_equal_zu ) THEN
-            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),50.0), zeta_t )
+            zeta_t = zt*ztmp0
+            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
          END IF
 
          !! Turbulent scales at zu=10m :
          ztmp0   = psi_h_coare(zeta_u)
-         ztmp1   = vkarmn/(LOG(zu) -LOG(z0t) - ztmp0)
+         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
 
          t_star = dt_zu*ztmp1
          q_star = dq_zu*ztmp1
-         u_star = U_blk*vkarmn/(LOG(zu) -LOG(z0) - psi_m_coare(zeta_u))
+         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
 
          IF( .NOT. l_zt_equal_zu ) THEN
@@ -259,5 +258,5 @@
       !! Wind greater than 18 m/s :  alfa = 0.018
       !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://sourceforge.net/p/aerobulk)
+      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
       !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
@@ -274,6 +273,6 @@
             !
             ! Charnock's constant, increases with the wind :
-            zgt10 = 0.5 + SIGN(0.5,(zw - 10.)) ! If zw<10. --> 0, else --> 1
-            zgt18 = 0.5 + SIGN(0.5,(zw - 18.)) ! If zw<18. --> 0, else --> 1
+            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1
+            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
             !
             alfa_charn(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
@@ -294,5 +293,5 @@
       !!
       !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://sourceforge.net/p/aerobulk)
+      !!         (https://github.com/brodeau/aerobulk/)
       !!------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
@@ -330,5 +329,5 @@
       !!       form from Beljaars and Holtslag (1991)
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
@@ -356,6 +355,6 @@
             zf = zta*zta
             zf = zf/(1. + zf)
-            zc = MIN(50., 0.35*zta)
-            zstab = 0.5 + SIGN(0.5, zta)
+            zc = MIN(50._wp, 0.35_wp*zta)
+            zstab = 0.5 + SIGN(0.5_wp, zta)
             !
             psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
@@ -383,5 +382,5 @@
       !!
       !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://sourceforge.net/p/aerobulk)
+      !!         (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------
       !!
@@ -408,6 +407,6 @@
             zf = zta*zta
             zf = zf/(1. + zf)
-            zc = MIN(50.,0.35*zta)
-            zstab = 0.5 + SIGN(0.5, zta)
+            zc = MIN(50._wp,0.35_wp*zta)
+            zstab = 0.5 + SIGN(0.5_wp, zta)
             !
             psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
@@ -420,29 +419,4 @@
    END FUNCTION psi_h_coare
 
-
-   FUNCTION visc_air( ptak )
-      !!---------------------------------------------------------------------
-      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
-      !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://sourceforge.net/p/aerobulk)
-      !!---------------------------------------------------------------------
-      !!
-      REAL(wp), DIMENSION(jpi,jpj) :: visc_air
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak  ! air temperature in (K)
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::   ztc, ztc2      ! local scalar
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
-            ztc2 = ztc*ztc
-            visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc)
-         END DO
-      END DO
-      !
-   END FUNCTION visc_air
-
-
+   !!======================================================================
 END MODULE sbcblk_algo_coare
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90	(revision 11110)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90	(revision 11111)
@@ -14,5 +14,5 @@
    !!
    !!       Routine turb_coare3p5 maintained and developed in AeroBulk
-   !!                     (http://aerobulk.sourceforge.net/)
+   !!                     (https://github.com/brodeau/aerobulk/)
    !!
    !!            Author: Laurent Brodeau, 2016, brodeau@gmail.com
@@ -51,5 +51,4 @@
    REAL(wp), PARAMETER ::   zi0        = 600.      ! scale height of the atmospheric boundary layer...1
    REAL(wp), PARAMETER ::   Beta0      =   1.25    ! gustiness parameter
-   REAL(wp), PARAMETER ::   rctv0      =   0.608   ! constant to obtain virtual temperature...
 
    !!----------------------------------------------------------------------
@@ -68,5 +67,5 @@
       !! ** Method : Monin Obukhov Similarity Theory
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
       !!
       !! INPUT :
@@ -265,5 +264,5 @@
       !!
       !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://sourceforge.net/p/aerobulk)
+      !!         (https://github.com/brodeau/aerobulk/)
       !!------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
@@ -300,5 +299,5 @@
       !!       form from Beljaars and Holtslag (1991)
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
@@ -353,5 +352,5 @@
       !!
       !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://sourceforge.net/p/aerobulk)
+      !!         (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------
       !!
@@ -396,5 +395,5 @@
       !!
       !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://sourceforge.net/p/aerobulk)
+      !!         (https://github.com/brodeau/aerobulk/)
       !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air
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 11110)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11111)
@@ -1,8 +1,6 @@
 MODULE sbcblk_algo_ecmwf
    !!======================================================================
-   !!                       ***  MODULE  sbcblk_algo_ecmwf  ***
-   !! Computes turbulent components of surface fluxes
-   !!         according to the method in IFS of the ECMWF model
-   !!
+   !!                   ***  MODULE  sbcblk_algo_ecmwf  ***
+   !! Computes:
    !!   * bulk transfer coefficients C_D, C_E and C_H
    !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
@@ -10,12 +8,11 @@
    !!   => all these are used in bulk formulas in sbcblk.F90
    !!
-   !!    Using the bulk formulation/param. of IFS of ECMWF (cycle 31r2)
+   !!    Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1)
    !!         based on IFS doc (avaible online on the ECMWF's website)
    !!
+   !!       Routine turb_ecmwf maintained and developed in AeroBulk
+   !!                     (https://github.com/brodeau/aerobulk)
    !!
-   !!       Routine turb_ecmwf maintained and developed in AeroBulk
-   !!                     (http://aerobulk.sourceforge.net/)
-   !!
-   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
@@ -41,4 +38,6 @@
 
    USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbcblk_phymbl   !LB: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF)
 
    IMPLICIT NONE
@@ -52,26 +51,48 @@
    REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1
    REAL(wp), PARAMETER ::   Beta0    = 1.     ! gustiness parameter ( = 1.25 in COAREv3)
-   REAL(wp), PARAMETER ::   rctv0    = 0.608  ! constant to obtain virtual temperature...
-   REAL(wp), PARAMETER ::   Cp_dry = 1005.0   ! Specic heat of dry air, constant pressure      [J/K/kg]
-   REAL(wp), PARAMETER ::   Cp_vap = 1860.0   ! Specic heat of water vapor, constant pressure  [J/K/kg]
    REAL(wp), PARAMETER ::   alpha_M = 0.11    ! For roughness length (smooth surface term)
    REAL(wp), PARAMETER ::   alpha_H = 0.40    ! (Chapter 3, p.34, IFS doc Cy31r1)
    REAL(wp), PARAMETER ::   alpha_Q = 0.62    !
+
+   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations
+
+   !! Cool-Skin / Warm-Layer related parameters:
+   REAL(wp),    PARAMETER :: &
+      &  rdt0    = 3600.*1.5, & !: time step
+      &  rd0     = 3. ,      &  !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005)
+      &  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 TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu,   &
-      &                   Cd, Ch, Ce , t_zu, q_zu, U_blk,         &
-      &                   Cdn, Chn, Cen                           )
+   SUBROUTINE TURB_ECMWF( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
+      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
+      &                   Cdn, Chn, Cen,                      &
+      &                   Qsw, rad_lw, slp                   )
       !!----------------------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_ecmwf  ***
       !!
-      !!            2015: L. Brodeau (brodeau@gmail.com)
-      !!
       !! ** Purpose :   Computes turbulent transfert coefficients of surface
-      !!                fluxes according to IFS doc. (cycle 31)
+      !!                fluxes according to IFS doc. (cycle 40)
       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
-      !!
-      !! ** Method : Monin Obukhov Similarity Theory
+      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
+      !!
+      !!                Applies the cool-skin warm-layer correction of the SST to T_s
+      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave
+      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
+      !!                are provided as (optional) arguments!
       !!
       !! INPUT :
@@ -80,9 +101,19 @@
       !!    *  zu   : height for wind speed (generally 10m)                   [m]
       !!    *  U_zu : scalar wind speed at 10m                                [m/s]
-      !!    *  sst  : SST                                                     [K]
       !!    *  t_zt : potential air temperature at zt                         [K]
-      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
       !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
       !!
+      !! INPUT/OUTPUT:
+      !! -------------
+      !!    *  T_s  : SST or skin temperature                                 [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)
+      !!              -> 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!):
+      !! ---------------
+      !!    *  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]
       !!
       !! OUTPUT :
@@ -93,15 +124,15 @@
       !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
       !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
-      !!    *  U_blk  : bulk wind at 10m                                      [m/s]
-      !!
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!    *  U_blk  : bulk wind speed at 10m                                [m/s]
+      !!
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
       REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                   [kg/kg]
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg]
+      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]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
@@ -113,17 +144,33 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
       !
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
+      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
+      !
       INTEGER :: j_itt
-      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
-      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
-      !
-      REAL(wp), DIMENSION(jpi,jpj) ::   u_star, t_star, q_star,   &
+      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::  &
+         &  u_star, t_star, q_star, &
          &  dt_zu, dq_zu,    &
          &  znu_a,           & !: Nu_air, Viscosity of air
          &  Linv,            & !: 1/L (inverse of Monin Obukhov length...
          &  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(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
       !
@@ -131,42 +178,50 @@
       IF( ABS(zu - zt) < 0.01 )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
 
+      !! Initialization for cool skin:
+      IF( l_use_skin ) THEN
+         zsst   = T_s    ! save the bulk SST
+         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
+      END IF
 
       !! First guess of temperature and humidity at height zu:
-      t_zu = MAX( t_zt , 0.0  )   ! who knows what's given on masked-continental regions...
-      q_zu = MAX( q_zt , 1.e-6)   !               "
+      t_zu = MAX( t_zt , 0.0_wp  )   ! who knows what's given on masked-continental regions...
+      q_zu = MAX( q_zt , 1.e-6_wp)   !               "
 
       !! Pot. temp. difference (and we don't want it to be 0!)
-      dt_zu = t_zu - sst   ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.e-6), dt_zu )
-      dq_zu = q_zu - ssq   ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.e-9), dq_zu )
+      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 )
 
       znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K)
 
-      ztmp2 = 0.5 * 0.5  ! initial guess for wind gustiness contribution
+      ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution
       U_blk = SQRT(U_zu*U_zu + ztmp2)
 
-      ! z0     = 0.0001
-      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0
+      ztmp2   = 10000._wp     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
       ztmp0   = LOG(zu*ztmp2)
       ztmp1   = LOG(10.*ztmp2)
-      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
-
-      z0     = charn0*u_star*u_star/grav + 0.11*znu_a/u_star
-      z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t !
-
-      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
-
-      ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd
+      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
+
+      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
+      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
+      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
+      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
+
+      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
 
       ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk )   ! Ribu = Bulk Richardson number
 
       !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2):
-      ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 )
+      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
       func_m = ztmp0*ztmp2 ! temporary array !!
-      !!             Ribu < 0                                 Ribu > 0   Beta = 1.25
-      func_h = (1.-ztmp1)*(func_m/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &  ! temporary array !!! func_h == zeta_u
-         &  +     ztmp1*(func_m*(1. + 27./9.*ztmp2/ztmp0))
+      func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  Ribu < 0 ! temporary array !!! func_h == zeta_u
+         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  Ribu > 0
+      !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
 
       !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
-      ztmp0   =        vkarmn/(LOG(zu*z0t) - psi_h_ecmwf(func_h))
+      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
 
       u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h))
@@ -176,16 +231,14 @@
       ! What's need to be done if zt /= zu:
       IF( .NOT. l_zt_equal_zu ) THEN
-         !
          !! First update of values at zu (or zt for wind)
          ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu)    ! zt*func_h/zu == zeta_t
-         ztmp1 = log(zt/zu) + ztmp0
+         ztmp1 = LOG(zt/zu) + ztmp0
          t_zu = t_zt - t_star/vkarmn*ztmp1
          q_zu = q_zt - q_star/vkarmn*ztmp1
-         q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity :
-
-         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
-         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
+         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
          !
-      ENDIF
+         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 )
+      END IF
 
 
@@ -195,16 +248,13 @@
       !! First guess of inverse of Monin-Obukov length (1/L) :
       ztmp0 = (1. + rctv0*q_zu)  ! the factor to apply to temp. to get virt. temp...
-      Linv  =  grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / ( u_star*u_star * t_zu*ztmp0 )
+      Linv  =  grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / MAX( u_star*u_star * t_zu*ztmp0 , 1.E-9 ) ! #LOLO
+      Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
       !! Functions such as  u* = U_blk*vkarmn/func_m
-      ztmp1 = zu + z0
-      ztmp0 = ztmp1*Linv
-      func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv)
-      func_h = LOG(ztmp1*z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(1./z0t*Linv)
-
+      ztmp0 = zu*Linv
+      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv)
+      func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
 
       !! ITERATION BLOCK
-      !! ***************
-
       DO j_itt = 1, nb_itt
 
@@ -213,9 +263,11 @@
 
          !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :
-         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3, p.33, IFS doc - Cy31r1
+         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
+         !! Note: it is slightly different that the L we would get with the usual
+         !! expression, as in coare algorithm or in 'mod_phymbl.f90' (One_on_L_MO())
+         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
          !! Update func_m with new Linv:
-         ztmp1 = zu + z0
-         func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp1*Linv) + psi_m_ecmwf(z0*Linv)
+         func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu!
 
          !! Need to update roughness lengthes:
@@ -223,13 +275,14 @@
          ztmp2  = u_star*u_star
          ztmp1  = znu_a/u_star
-         z0    = alpha_M*ztmp1 + charn0*ztmp2/grav
-         z0t    = alpha_H*ztmp1                              ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
-         z0q    = alpha_Q*ztmp1
-
+         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
+         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
+         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.))**(2./3.) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
+         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)              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1
+         U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp )              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1
          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
 
@@ -238,12 +291,11 @@
          !! as well the air-sea differences:
          IF( .NOT. l_zt_equal_zu ) THEN
-
             !! Arrays func_m and func_h are free for a while so using them as temporary arrays...
-            func_h = psi_h_ecmwf((zu+z0)*Linv) ! temporary array !!!
-            func_m = psi_h_ecmwf((zt+z0)*Linv) ! temporary array !!!
+            func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!!
+            func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!!
 
             ztmp2  = psi_h_ecmwf(z0t*Linv)
             ztmp0  = func_h - ztmp2
-            ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0t) - ztmp0)
+            ztmp1  = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
             t_star = dt_zu*ztmp1
             ztmp2  = ztmp0 - func_m + ztmp2
@@ -253,32 +305,49 @@
             ztmp2  = psi_h_ecmwf(z0q*Linv)
             ztmp0  = func_h - ztmp2
-            ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0q) - ztmp0)
+            ztmp1  = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0)
             q_star = dq_zu*ztmp1
             ztmp2  = ztmp0 - func_m + ztmp2
-            ztmp1  = log(zt/zu) + ztmp2
+            ztmp1  = LOG(zt/zu) + ztmp2
             q_zu   = q_zt - q_star/vkarmn*ztmp1
-
-            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
-            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
-
          END IF
 
          !! Updating because of updated z0 and z0t and new Linv...
-         ztmp1 = zu + z0
-         ztmp0 = ztmp1*Linv
-         func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
-         func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
-
-      END DO
+         ztmp0 = zu*Linv
+         func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
+         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...
+            Cd = vkarmn*vkarmn/(func_m*func_m)
+            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*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
+            !! 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
+            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 )
+         END IF
+
+      END DO !DO j_itt = 1, nb_itt
 
       Cd = vkarmn*vkarmn/(func_m*func_m)
       Ch = vkarmn*vkarmn/(func_m*func_h)
-      ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
-      Ce = vkarmn*vkarmn/(func_m*ztmp1)
-
-      ztmp1 = zu + z0
-      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
-      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
-      Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q))
+      ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
+      Ce = vkarmn*vkarmn/(func_m*ztmp2)
+
+      Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 ))
+      Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t))
+      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
 
    END SUBROUTINE TURB_ECMWF
@@ -294,5 +363,5 @@
       !!         and L is M-O length
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf
@@ -306,24 +375,24 @@
          DO ji = 1, jpi
             !
-            zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!):
+            zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
             !
             ! Unstable (Paulson 1970):
             !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
-            zx = SQRT(ABS(1. - 16.*zzeta))
-            ztmp = 1. + SQRT(zx)
+            zx = SQRT(ABS(1._wp - 16._wp*zzeta))
+            ztmp = 1._wp + SQRT(zx)
             ztmp = ztmp*ztmp
-            psi_unst = LOG( 0.125*ztmp*(1. + zx) )   &
-               &       -2.*ATAN( SQRT(zx) ) + 0.5*rpi
+            psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
+               &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
             !
             ! Unstable:
             ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
-            psi_stab = -2./3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) &
-               &       - zzeta - 2./3.*5./0.35
+            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
+               &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
             !
             ! Combining:
-            stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1
-            !
-            psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable
-               &                +      stab  * psi_stab   ! (zzeta > 0) Stable
+            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
+            !
+            psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
+               &                +      stab  * psi_stab      ! (zzeta > 0) Stable
             !
          END DO
@@ -332,5 +401,5 @@
    END FUNCTION psi_m_ecmwf
 
-   
+
    FUNCTION psi_h_ecmwf( pzeta )
       !!----------------------------------------------------------------------------------
@@ -342,5 +411,5 @@
       !!         and L is M-O length
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
@@ -354,21 +423,21 @@
          DO ji = 1, jpi
             !
-            zzeta = MIN(pzeta(ji,jj) , 5.)   ! Very stable conditions (L positif and big!):
-            !
-            zx  = ABS(1. - 16.*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
+            zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
+            !
+            zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
             !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
             ! Unstable (Paulson 1970) :
-            psi_unst = 2.*LOG(0.5*(1. + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
+            psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
             !
             ! Stable:
-            psi_stab = -2./3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
-               &       - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1. 
+            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
+               &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
             ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
             !
-            stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1
-            !
-            !
-            psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst &   ! (zzeta < 0) Unstable
-               &                +    stab    * psi_stab     ! (zzeta > 0) Stable
+            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
+            !
+            !
+            psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
+               &                +    stab    * psi_stab        ! (zzeta > 0) Stable
             !
          END DO
@@ -382,5 +451,5 @@
       !! Bulk Richardson number (Eq. 3.25 IFS doc)
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj) ::   Ri_bulk   !
@@ -394,6 +463,6 @@
       !!----------------------------------------------------------------------------------
       !
-      Ri_bulk =   grav*pz/(pub*pub)                                          &
-         &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(Cp_dry+Cp_vap*pqz)))   &
+      Ri_bulk =   grav*pz/(pub*pub)   &
+         &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(rCp_dry + rCp_vap*pqz))) &
          &          + rctv0*pdq )
       !
@@ -401,26 +470,4 @@
 
 
-   FUNCTION visc_air(ptak)
-      !!----------------------------------------------------------------------------------
-      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
-      !
-      INTEGER  ::   ji, jj      ! dummy loop indices
-      REAL(wp) ::   ztc, ztc2   ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
-            ztc2 = ztc*ztc
-            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
-         END DO
-      END DO
-      !
-   END FUNCTION visc_air
 
    !!======================================================================
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11110)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11111)
@@ -11,5 +11,5 @@
    !!
    !!       Routine turb_ncar maintained and developed in AeroBulk
-   !!                     (http://aerobulk.sourceforge.net/)
+   !!                     (https://github.com/brodeau/aerobulk/)
    !!
    !!                         L. Brodeau, 2015
@@ -43,7 +43,4 @@
 
    PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90
-
-   !                              ! NCAR own values for given constants:
-   REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature...
    
    !!----------------------------------------------------------------------
@@ -76,5 +73,5 @@
       !!      => 'vkarmn' and 'grav'
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!
       !! INPUT :
@@ -238,5 +235,5 @@
       !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
@@ -277,5 +274,5 @@
       !!         and L is M-O length
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pzeta
@@ -311,5 +308,5 @@
       !!         and L is M-O length
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phymbl.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phymbl.F90	(revision 11111)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phymbl.F90	(revision 11111)
@@ -0,0 +1,193 @@
+MODULE sbcblk_phymbl
+   !!======================================================================
+   !!                       ***  MODULE  sbcblk_phymbl  ***
+   !! A set of functions to compute air themodynamics parameters
+   !!                     needed by Aerodynamic Bulk Formulas
+   !!=====================================================================
+   !!            4.0  !  2019 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
+   !!----------------------------------------------------------------------
+
+   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
+   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
+   !!   q_sat         : saturation humidity as a function of SLP and temperature
+   !!   gamma_moist   : 
+   !!   L_vap         : latent heat of vaporization of water as a function of temperature
+   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature
+   
+   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 sbcdcy         ! surface boundary condition: diurnal cycle
+   !USE sbcwave , ONLY :   cdn_wave ! wave module
+   !USE sbc_ice        ! Surface boundary condition: ice fields
+   USE lib_fortran    ! to use key_nosignedzero
+   !
+   !USE iom            ! I/O manager library
+   !USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing library
+   !USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   !USE prtctl         ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC rho_air
+   PUBLIC cp_air
+   PUBLIC q_sat
+   PUBLIC gamma_moist
+   PUBLIC L_vap
+   PUBLIC visc_air
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   FUNCTION rho_air( ptak, pqa, pslp )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION rho_air  ***
+      !!
+      !! ** 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), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
+      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
+
+
+   FUNCTION cp_air( pqa )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air  ***
+      !!
+      !! ** Purpose : Compute specific heat (Cp) of moist air
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!-------------------------------------------------------------------------------
+      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
+
+
+   FUNCTION q_sat( ptak, pslp )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION q_sat  ***
+      !!
+      !! ** Purpose : Specific humidity at saturation in [kg/kg]
+      !!              Based on accurate estimate of "e_sat"
+      !!              aka saturation water vapor (Goff, 1957)
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) ::   ze_sat, ztmp   ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            ztmp = rt0 / ptak(ji,jj)
+            !
+            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
+            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
+               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
+               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
+            !
+            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION q_sat
+
+
+   FUNCTION gamma_moist( ptak, pqa )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION gamma_moist  ***
+      !!
+      !! ** Purpose : Compute the moist adiabatic lapse-rate.
+      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
+      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) :: zrv, ziRT        ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
+            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
+            gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( rCp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )
+         END DO
+      END DO
+      !
+   END FUNCTION gamma_moist
+
+
+   FUNCTION L_vap( psst )
+      !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION L_vap  ***
+      !!
+      !! ** 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   ! 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 visc_air(ptak)
+      !!----------------------------------------------------------------------------------
+      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
+      !
+      INTEGER  ::   ji, jj      ! dummy loop indices
+      REAL(wp) ::   ztc, ztc2   ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
+            ztc2 = ztc*ztc
+            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
+         END DO
+      END DO
+      !
+   END FUNCTION visc_air
+   
+
+   !!======================================================================
+END MODULE sbcblk_phymbl
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 11111)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin.F90	(revision 11111)
@@ -0,0 +1,223 @@
+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_phymbl  !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_*.F90
+
+   !! Cool-Skin / Warm-Layer related parameters:
+   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   = rdt0/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
