Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11139)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 12928)
@@ -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 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
@@ -41,48 +38,109 @@
 
    USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
+   USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC ::   TURB_ECMWF   ! called by sbcblk.F90
-
-   !                   !! ECMWF own values for given constants, taken form IFS documentation...
+   PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+
+   !! ECMWF own values for given constants, taken form IFS documentation...
    REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!!
    !                                          !    =>  Usually 0.011 for moderate winds)
    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 = 10             ! number of itterations
+
    !!----------------------------------------------------------------------
 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                           )
-      !!----------------------------------------------------------------------------------
-      !!                      ***  ROUTINE  turb_ecmwf  ***
-      !!
-      !!            2015: L. Brodeau (brodeau@gmail.com)
-      !!
-      !! ** Purpose :   Computes turbulent transfert coefficients of surface
-      !!                fluxes according to IFS doc. (cycle 31)
-      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
-      !!
-      !! ** Method : Monin Obukhov Similarity Theory
+
+   SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl)
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION sbcblk_algo_ecmwf_init  ***
       !!
       !! INPUT :
       !! -------
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
+      !!---------------------------------------------------------------------
+      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization
+      INTEGER :: ierr
+      !!---------------------------------------------------------------------
+      IF( l_use_wl ) THEN
+         ierr = 0
+         ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr )
+         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' )
+         dT_wl(:,:)  = 0._wp
+         Hz_wl(:,:)  = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars)
+      ENDIF
+      IF( l_use_cs ) THEN
+         ierr = 0
+         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
+         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' )
+         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
+      ENDIF
+   END SUBROUTINE sbcblk_algo_ecmwf_init
+
+
+
+   SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
+      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                           &
+      &                      Cdn, Chn, Cen,                                           &
+      &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer)
+      &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only
+      !!----------------------------------------------------------------------------------
+      !!                      ***  ROUTINE  turb_ecmwf  ***
+      !!
+      !! ** Purpose :   Computes turbulent transfert coefficients of surface
+      !!                fluxes according to IFS doc. (cycle 45r1)
+      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
+      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
+      !!
+      !!                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 :
+      !! -------
+      !!    *  kt   : current time step (starts at 1)
       !!    *  zt   : height for temperature and spec. hum. of air            [m]
-      !!    *  zu   : height for wind speed (generally 10m)                   [m]
-      !!    *  U_zu : scalar wind speed at 10m                                [m/s]
-      !!    *  sst  : SST                                                     [K]
+      !!    *  zu   : height for wind speed (usually 10m)                     [m]
       !!    *  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]
-      !!
+      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
+      !!
+      !! INPUT/OUTPUT:
+      !! -------------
+      !!    *  T_s  : always "bulk SST" as input                              [K]
+      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
+      !!              -> skin temperature as output if CSWL used              [K]
+      !!
+      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
+      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True)
+      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False)
+      !!
+      !! OPTIONAL INPUT:
+      !! ---------------
+      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
+      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
+      !!    *  slp    : sea-level pressure                                    [Pa]
+      !!
+      !! OPTIONAL OUTPUT:
+      !! ----------------
+      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
+      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
+      !!    * pHz_wl  : thickness of warm-layer                               [m]
       !!
       !! OUTPUT :
@@ -93,16 +151,19 @@
       !!    *  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 zu                                 [m/s]
+      !!
+      !!
+      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
       REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
-      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]
+      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
@@ -110,81 +171,103 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s]
       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]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m]
+      !
       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,   &
-         &  dt_zu, dq_zu,    &
-         &  znu_a,           & !: Nu_air, Viscosity of air
-         &  Linv,            & !: 1/L (inverse of Monin Obukhov length...
-         &  z0, z0t, z0q
-      REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
-      !!----------------------------------------------------------------------------------
-      !
-      ! Identical first gess as in COARE, with IFS parameter values though
-      !
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01 )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
-
-
+      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
+      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu
+      REAL(wp), DIMENSION(jpi,jpj) :: znu_a         !: Nu_air, Viscosity of air
+      REAL(wp), DIMENSION(jpi,jpj) :: Linv  !: 1/L (inverse of Monin Obukhov length...
+      REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q
+      !
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: 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
+      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
+      !!----------------------------------------------------------------------------------
+      IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl)
+
+      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
+
+      !! Initializations for cool skin and warm layer:
+      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
+
+      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
+
+      IF( l_use_cs .OR. l_use_wl ) THEN
+         ALLOCATE ( zsst(jpi,jpj) )
+         zsst = T_s ! backing up the bulk SST
+         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
+      ENDIF
+
+
+      ! Identical first gess as in COARE, with IFS parameter values though...
+      !
       !! First guess of temperature and humidity at height zu:
-      t_zu = MAX( t_zt , 0.0  )   ! who knows what's given on masked-continental regions...
-      q_zu = MAX( q_zt , 1.e-6)   !               "
+      t_zu = MAX( t_zt ,  180._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)
-
-      ztmp2 = 0.5 * 0.5  ! initial guess for wind gustiness contribution
-      U_blk = SQRT(U_zu*U_zu + ztmp2)
-
-      ! z0     = 0.0001
-      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0
-      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 !
+      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_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
+
+      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
+
+      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
+      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
+      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( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
+
+      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
+      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
 
       Cd     = (vkarmn/ztmp0)**2    ! 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 )
+      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
+
+      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
+
+      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (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._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u
+         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 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))
-
-      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h))
+      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
+
+      u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on)
       t_star = dt_zu*ztmp0
       q_star = dq_zu*ztmp0
 
-      ! What's need to be done if zt /= zu:
+      ! What needs to be done if zt /= zu:
       IF( .NOT. l_zt_equal_zu ) THEN
-         !
          !! First update of values at zu (or zt for wind)
          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 :
          !
+         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 )
       ENDIF
 
@@ -194,28 +277,24 @@
 
       !! 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 = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
 
       !! 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
 
          !! Bulk Richardson Number at z=zu (Eq. 3.25)
-         ztmp0 = Ri_bulk(zu, t_zu, dt_zu, q_zu, dq_zu, U_blk)
+         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
 
          !! 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
+         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...)
 
          !! 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 +302,12 @@
          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
-
-         !! Update wind at 10m taking into acount convection-related wind gustiness:
-         ! 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)
-         !! => 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
+         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 zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8)
+         ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
+         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
+         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
 
@@ -238,12 +316,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,34 +330,65 @@
             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
+         ENDIF
 
          !! 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)
+
+
+         IF( l_use_cs ) THEN
+            !! Cool-skin contribution
+
+            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, &
+               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0
+
+            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst )  ! Qnsol -> ztmp1
+
+            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1)
+            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+
+         ENDIF
+
+         IF( l_use_wl ) THEN
+            !! Warm-layer contribution
+            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, &
+               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2
+            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst )
+            !! Updating T_s and q_s !!!
+            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) !
+            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+         ENDIF
+
+         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
+            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
+            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
+         ENDIF
+
+      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))
-
-   END SUBROUTINE TURB_ECMWF
+      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))
+
+      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
+      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
+      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
+
+      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+
+   END SUBROUTINE turb_ecmwf
 
 
@@ -294,5 +402,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
@@ -302,35 +410,31 @@
       REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab
       !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zzeta = MIN( pzeta(ji,jj) , 5. ) !! 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)
-            ztmp = ztmp*ztmp
-            psi_unst = LOG( 0.125*ztmp*(1. + zx) )   &
-               &       -2.*ATAN( SQRT(zx) ) + 0.5*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
-            !
-            ! 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
-            !
-         END DO
-      END DO
-      !
+      DO_2D_11_11
+      !
+      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._wp - 16._wp*zzeta))
+      ztmp = 1._wp + SQRT(zx)
+      ztmp = ztmp*ztmp
+      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._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_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_2D
    END FUNCTION psi_m_ecmwf
 
-   
+
    FUNCTION psi_h_ecmwf( pzeta )
       !!----------------------------------------------------------------------------------
@@ -342,5 +446,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
@@ -351,76 +455,27 @@
       !!----------------------------------------------------------------------------------
       !
-      DO jj = 1, jpj
-         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  !!!
-            !                                     ! 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
-            !
-            ! 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. 
-            ! 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
-            !
-         END DO
-      END DO
-      !
+      DO_2D_11_11
+      !
+      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._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
+      !
+      ! Stable:
+      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_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_2D
    END FUNCTION psi_h_ecmwf
 
-
-   FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub )
-      !!----------------------------------------------------------------------------------
-      !! Bulk Richardson number (Eq. 3.25 IFS doc)
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) ::   Ri_bulk   !
-      !
-      REAL(wp)                    , INTENT(in) ::   pz    ! height above the sea        [m]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptz   ! air temperature at pz m     [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdt   ! ptz - sst                   [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqz   ! air temperature at pz m [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdq   ! pqz - ssq               [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pub   ! bulk wind speed           [m/s]
-      !!----------------------------------------------------------------------------------
-      !
-      Ri_bulk =   grav*pz/(pub*pub)                                          &
-         &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(Cp_dry+Cp_vap*pqz)))   &
-         &          + rctv0*pdq )
-      !
-   END FUNCTION Ri_bulk
-
-
-   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
 
    !!======================================================================
