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 11283)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11284)
@@ -46,6 +46,6 @@
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
-   USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003)
-   USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
+   USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
+   USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013)
    USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)
    !
@@ -87,5 +87,5 @@
    LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
    LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
-   LOGICAL  ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013)
+   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013)
    LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
    !
@@ -109,5 +109,13 @@
    !LB:
    LOGICAL  ::   ln_skin        ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
-   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature if .true. (specific humidity espected if .false.) !LB
+   LOGICAL  ::   ln_humi_sph    ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB
+   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
+   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB
+   !
+   INTEGER  ::   nhumi          ! choice of the bulk algorithm
+   !                            ! associated indices:
+   INTEGER, PARAMETER :: np_humi_sph = 1
+   INTEGER, PARAMETER :: np_humi_dpt = 2
+   INTEGER, PARAMETER :: np_humi_rlh = 3
    !LB.
 
@@ -116,5 +124,5 @@
    INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008)
    INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
-   INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013)
+   INTEGER, PARAMETER ::   np_COARE_3p6 = 3   ! "COARE 3.6" algorithm   (Edson et al. 2013)
    INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31)
 
@@ -172,8 +180,8 @@
       NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
-         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
+         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm
          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &
          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
-         &                 ln_skin, ln_humi_dpt                                                        ! cool-skin / warm-layer param. !LB
+         &                 ln_skin, ln_humi_sph, ln_humi_dpt, ln_humi_rlh    ! cool-skin / warm-layer !LB
       !!---------------------------------------------------------------------
       !
@@ -201,11 +209,14 @@
          nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1
       ENDIF
-      IF( ln_COARE_3p5 ) THEN
-         nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1
+      IF( ln_COARE_3p6 ) THEN
+         nblk =  np_COARE_3p6   ;   ioptio = ioptio + 1
       ENDIF
       IF( ln_ECMWF     ) THEN
          nblk =  np_ECMWF       ;   ioptio = ioptio + 1
       ENDIF
-      !
+      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
+
+
+
       !LB:
       !                             !** initialization of the cool-skin / warm-layer parametrization
@@ -215,7 +226,18 @@
          IF( sbc_blk_cswl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
       END IF
+      !
+      ioptio = 0
+      IF( ln_humi_sph ) THEN
+         nhumi =  np_humi_sph    ;   ioptio = ioptio + 1
+      ENDIF
+      IF( ln_humi_dpt ) THEN
+         nhumi =  np_humi_dpt    ;   ioptio = ioptio + 1
+      ENDIF
+      IF( ln_humi_rlh ) THEN
+         nhumi =  np_humi_rlh    ;   ioptio = ioptio + 1
+      ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' )
       !LB.
 
-      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
       !
       IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
@@ -278,5 +300,5 @@
          WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
          WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
-         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p5
+         WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6
          WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
          WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
@@ -294,5 +316,5 @@
          CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)'
          CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
-         CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)'
+         CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
          CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)'
          END SELECT
@@ -300,6 +322,13 @@
          WRITE(numout,*)
          WRITE(numout,*) '      use cool-skin/warm-layer parameterization (SSST)     ln_skin     = ', ln_skin !LB
+         !
+         !LB:
          WRITE(numout,*)
-         WRITE(numout,*) '      air humidity read in files is dew-point temperature? ln_humi_dpt = ', ln_humi_dpt !LB
+         SELECT CASE( nhumi )              !* Print the choice of air humidity
+         CASE( np_humi_sph )   ;   WRITE(numout,*) '   ==>>>   air humidity is SPECIFIC HUMIDITY     [kg/kg]'
+         CASE( np_humi_dpt )   ;   WRITE(numout,*) '   ==>>>   air humidity is DEW-POINT TEMPERATURE [K]'
+         CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]'
+         END SELECT
+         !LB.
          !
       ENDIF
@@ -360,9 +389,12 @@
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
          !LB:
-         IF ( ln_humi_dpt ) THEN
-            qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
-         ELSE
-            qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1)
-         END IF
+         SELECT CASE( nhumi )
+         CASE( np_humi_sph )
+            qatm_ice(:,:) =           sf(jp_humi)%fnow(:,:,1)
+         CASE( np_humi_dpt )
+            qatm_ice(:,:) = q_sat(    sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+         CASE( np_humi_rlh )
+            qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !RH is % percent in file
+         END SELECT
          !LB.
          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
@@ -471,12 +503,15 @@
 
       !LB:
-      ! spec. hum. of air at "rn_zqt" m above thes sea:
-      IF ( ln_humi_dpt ) THEN
-         ! spec. hum. of air is deduced from dew-point temperature and SLP: q_air = q_sat( d_air, slp)
+      ! zqair = specific humidity of air at "rn_zqt" m above thes sea:
+      SELECT CASE( nhumi )
+      CASE( np_humi_sph )
+         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
+      CASE( np_humi_dpt )
          IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !'
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
-      ELSE
-         zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity!
-      END IF
+      CASE( np_humi_rlh )
+         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of RH, t_air and slp !'
+         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+      END SELECT
       !LB.
 
@@ -485,12 +520,18 @@
       !!    (since reanalysis products provide T at z, not theta !)
       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt
-      
+
 
       IF ( ln_skin ) THEN
-         
+
          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
 
          CASE( np_COARE_3p0 )
-            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
+            CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
+               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
+               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
+               &                Tsk_b=tsk(:,:) )
+
+         CASE( np_COARE_3p6 )
+            CALL turb_coare3p6 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, &
                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
@@ -502,5 +543,5 @@
                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
                &                Tsk_b=tsk(:,:) )
-            
+
          CASE DEFAULT
             CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' )
@@ -523,5 +564,5 @@
       ELSE !IF ( ln_skin )
 
-         
+
          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
             !
@@ -530,9 +571,9 @@
 
          CASE( np_COARE_3p0 )
-            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!'
-            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
+            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!'
+            CALL turb_coare3p0   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
-         CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5
+         CASE( np_COARE_3p6 )   ;   CALL turb_coare3p6( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6
             &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
@@ -802,12 +843,14 @@
 
       !LB:
-      ! spec. hum. of air at "rn_zqt" m above thes sea:
-      IF ( ln_humi_dpt ) THEN
-         ! spec. hum. of air is deduced from dew-point temperature and SLP: q_air = q_sat( d_air, slp)
+      SELECT CASE( nhumi )
+      CASE( np_humi_sph )
+         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
+      CASE( np_humi_dpt )
          IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !'
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
-      ELSE
-         zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity!
-      END IF
+      CASE( np_humi_rlh )
+         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !'
+         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+      END SELECT
       !LB.
 
Index: MO/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 11283)
+++ 	(revision )
@@ -1,457 +1,0 @@
-MODULE sbcblk_algo_coare
-   !!======================================================================
-   !!                   ***  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
-   !!   * the effective bulk wind speed at 10m U_blk
-   !!   => all these are used in bulk formulas in sbcblk.F90
-   !!
-   !!    Using the bulk formulation/param. of COARE v3, Fairall et al., 2003
-   !!
-   !!       Routine turb_coare 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
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   turb_coare  : computes the bulk turbulent transfer coefficients
-   !!                   adjusts t_air and q_air from zt to zu m
-   !!                   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 iom             ! I/O manager library
-   USE lib_mpp         ! distribued memory computing library
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE sbcwave, ONLY   :  cdn_wave ! wave module
-#if defined key_si3 || defined key_cice
-   USE sbc_ice         ! Surface boundary condition: ice fields
-#endif
-   USE lib_fortran     ! to use key_nosignedzero
-
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
-   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
-
-   !                                              !! COARE own values for given constants:
-   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
-   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
-
-   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
-
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE turb_coare( 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, Tsk_b             )
-      !!----------------------------------------------------------------------
-      !!                      ***  ROUTINE  turb_coare  ***
-      !!
-      !! ** 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
-      !!                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 :
-      !! -------
-      !!    *  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]
-      !!    *  t_zt : potential air temperature at zt                         [K]
-      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
-      !!
-      !! 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_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]
-      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
-      !!
-      !! OUTPUT :
-      !! --------
-      !!    *  Cd     : drag coefficient
-      !!    *  Ch     : sensible heat coefficient
-      !!    *  Ce     : evaporation coefficient
-      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
-      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
-      !!    *  U_blk  : bulk wind speed at 10m                                [m/s]
-      !!
-      !!
-      !! ** Author: L. Brodeau, June 2019 / 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(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(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)
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
-      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) ::   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(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
-      !
-      INTEGER :: j_itt
-      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
-         &  z0, z0t
-      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
-      !
-      ! Cool skin:
-      LOGICAL :: l_use_skin = .FALSE.
-      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
-      !!----------------------------------------------------------------------------------
-      !
-      ! 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
-
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
-
-      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
-
-      !! Initialization for cool skin:
-      zsst   = T_s    ! save the bulk SST
-      IF( l_use_skin ) THEN
-         ! First guess for skin temperature:
-         IF( PRESENT(Tsk_b) ) THEN
-            T_s = Tsk_b
-         ELSE
-            T_s = T_s - 0.25     ! sst - 0.25
-         END IF
-         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
-      END IF
-
-      !! First guess of temperature and humidity at height zu:
-      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
-      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 - 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)
-
-      ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution
-      U_blk = SQRT(U_zu*U_zu + ztmp2)
-
-      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_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
-
-      z0     = alfa_charn(U_blk)*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_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 )
-      ztmp0 = ztmp0*ztmp2
-      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
-         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
-      !#LB: 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))
-
-      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
-      t_star = dt_zu*ztmp0
-      q_star = dq_zu*ztmp0
-
-      ! 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)
-         zeta_t = zt*zeta_u/zu
-         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
-         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_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 )
-      END IF
-
-      !! ITERATION BLOCK
-      DO j_itt = 1, nb_itt
-
-         !!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
-
-         !! 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._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_wp)        ! include gustiness in bulk wind speed
-         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
-
-         !! Updating Charnock parameter, increases with the wind (Fairall et al., 2003 p. 577-578)
-         ztmp2 = alfa_charn(U_blk)  ! alpha Charnock parameter
-
-         !! Roughness lengthes z0, z0t (z0q = z0t) :
-         z0    = ztmp2*ztmp1/grav + 0.11_wp*znu_a/u_star  ! Roughness length (eq.6)
-         ztmp1 = z0*u_star/znu_a                          ! Re_r: roughness Reynolds number
-         z0t  = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28)
-
-         !! Stability parameters:
-         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_wp), zeta_t )
-         END IF
-
-         !! Turbulent scales at zu=10m :
-         ztmp0   = psi_h_coare(zeta_u)
-         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))
-
-         IF( .NOT. l_zt_equal_zu ) THEN
-            ! What's need to be done if zt /= zu
-            !! Re-updating temperature and humidity at zu :
-            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
-            ztmp1 = log(zt/zu) + ztmp2
-            t_zu = t_zt - t_star/vkarmn*ztmp1
-            q_zu = q_zt - q_star/vkarmn*ztmp1
-         END IF
-
-         !! SKIN related part
-         !! -----------------
-         IF( l_use_skin ) THEN
-            !! compute transfer coefficients at zu : lolo: verifier...
-            ztmp0 = u_star/U_blk
-            Ch   = ztmp0*t_star/dt_zu
-            Ce   = ztmp0*q_star/dq_zu
-            ! Non-Solar heat flux to the ocean:
-            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
-            ztmp2 = T_s*T_s
-            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
-               &       +      rad_lw - emiss_w*stefan*ztmp2*ztmp2                           ! Net longwave flux
-            !!         => "ztmp1" is the net non-solar surface heat flux !
-            !! Updating the values of the skin temperature T_s and q_s :
-            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
-            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
-
-      ! compute transfer coefficients at zu :
-      ztmp0 = u_star/U_blk
-      Cd   = ztmp0*ztmp0
-      Ch   = ztmp0*t_star/dt_zu
-      Ce   = ztmp0*q_star/dq_zu
-      !
-      ztmp1 = zu + z0
-      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
-      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
-      Cen = Chn
-      !
-      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
-      !
-   END SUBROUTINE turb_coare
-
-
-   FUNCTION alfa_charn( pwnd )
-      !!-------------------------------------------------------------------
-      !! Compute the Charnock parameter as a function of the wind speed
-      !!
-      !! (Fairall et al., 2003 p.577-578)
-      !!
-      !! Wind below 10 m/s :  alfa = 0.011
-      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
-      !! Wind greater than 18 m/s :  alfa = 0.018
-      !!
-      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
-      !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) :: zw, zgt10, zgt18
-      !!-------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zw = pwnd(ji,jj)   ! wind speed
-            !
-            ! Charnock's constant, increases with the wind :
-            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
-               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
-               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION alfa_charn
-
-   FUNCTION psi_m_coare( pzeta )
-      !!----------------------------------------------------------------------------------
-      !! ** Purpose: compute the universal profile stability function for momentum
-      !!             COARE 3.0, Fairall et al. 2003
-      !!             pzeta : stability paramenter, z/L where z is altitude
-      !!                     measurement and L is M-O length
-      !!       Stability function for wind speed and scalars matching Kansas and free
-      !!       convection forms with weighting f convective form, follows Fairall et
-      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
-      !!       form from Beljaars and Holtslag (1991)
-      !!
-      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zta = pzeta(ji,jj)
-            !
-            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
-            !
-            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
-               & - 2.*ATAN(zphi_m) + 0.5*rpi
-            !
-            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
-            !
-            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-            !
-            zf = zta*zta
-            zf = zf/(1. + zf)
-            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)
-               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
-               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION psi_m_coare
-
-
-   FUNCTION psi_h_coare( pzeta )
-      !!---------------------------------------------------------------------
-      !! Universal profile stability function for temperature and humidity
-      !! COARE 3.0, Fairall et al. 2003
-      !!
-      !! pzeta : stability paramenter, z/L where z is altitude measurement
-      !!         and L is M-O length
-      !!
-      !! Stability function for wind speed and scalars matching Kansas and free
-      !! convection forms with weighting f convective form, follows Fairall et
-      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
-      !! form from Beljaars and Holtslag (1991)
-      !!
-      !! Author: L. Brodeau, June 2016 / AeroBulk
-      !!         (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------
-      !!
-      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
-      !
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zta = pzeta(ji,jj)
-            !
-            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
-            !
-            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
-            !
-            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
-            !
-            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-            !
-            zf = zta*zta
-            zf = zf/(1. + zf)
-            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 ) &
-               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
-               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION psi_h_coare
-
-   !!======================================================================
-END MODULE sbcblk_algo_coare
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11284)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11284)
@@ -0,0 +1,456 @@
+MODULE sbcblk_algo_coare3p0
+   !!======================================================================
+   !!                   ***  MODULE  sbcblk_algo_coare3p0  ***
+   !! 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
+   !!   * the effective bulk wind speed at 10m U_blk
+   !!   => all these are used in bulk formulas in sbcblk.F90
+   !!
+   !!    Using the bulk formulation/param. of COARE v3, Fairall et al., 2003
+   !!
+   !!       Routine turb_coare3p0 maintained and developed in AeroBulk
+   !!                     (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
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   turb_coare3p0  : computes the bulk turbulent transfer coefficients
+   !!                   adjusts t_air and q_air from zt to zu m
+   !!                   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 iom             ! I/O manager library
+   USE lib_mpp         ! distribued memory computing library
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+   USE sbcwave, ONLY   :  cdn_wave ! wave module
+#if defined key_si3 || defined key_cice
+   USE sbc_ice         ! Surface boundary condition: ice fields
+#endif
+   USE lib_fortran     ! to use key_nosignedzero
+
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
+   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC ::   TURB_COARE3P0   ! called by sbcblk.F90
+
+   !                                              !! COARE own values for given constants:
+   REAL(wp), PARAMETER ::   zi0     = 600._wp     ! scale height of the atmospheric boundary layer...
+   REAL(wp), PARAMETER ::   Beta0   =   1.25_wp   ! gustiness parameter
+
+   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
+
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE turb_coare3p0( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
+      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
+      &                      Cdn, Chn, Cen,                      &
+      &                      Qsw, rad_lw, slp, Tsk_b             )
+      !!----------------------------------------------------------------------
+      !!                      ***  ROUTINE  turb_coare3p0  ***
+      !!
+      !! ** 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
+      !!                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 :
+      !! -------
+      !!    *  zt   : height for temperature and spec. hum. of air            [m]
+      !!    *  zu   : height for wind speed (usually 10m)                     [m]
+      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
+      !!    *  t_zt : potential air temperature at zt                         [K]
+      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
+      !!
+      !! 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_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]
+      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
+      !!
+      !! OUTPUT :
+      !! --------
+      !!    *  Cd     : drag coefficient
+      !!    *  Ch     : sensible heat coefficient
+      !!    *  Ce     : evaporation coefficient
+      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
+      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
+      !!    *  U_blk  : bulk wind speed at zu                                 [m/s]
+      !!
+      !!
+      !! ** Author: L. Brodeau, June 2019 / 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(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(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)
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
+      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 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(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
+      !
+      INTEGER :: j_itt
+      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
+         &  z0, z0t
+      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
+      !
+      ! Cool skin:
+      LOGICAL :: l_use_skin = .FALSE.
+      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
+      !!----------------------------------------------------------------------------------
+      !
+      ! 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
+
+      l_zt_equal_zu = .FALSE.
+      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+
+      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
+
+      !! Initialization for cool skin:
+      zsst   = T_s    ! save the bulk SST
+      IF( l_use_skin ) THEN
+         ! First guess for skin temperature:
+         IF( PRESENT(Tsk_b) ) THEN
+            T_s = Tsk_b
+         ELSE
+            T_s = T_s - 0.25     ! sst - 0.25
+         END IF
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
+      END IF
+
+      !! First guess of temperature and humidity at height zu:
+      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
+      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 - 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     = alfa_charn_3p0(U_zu)*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_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 )
+      ztmp0 = ztmp0*ztmp2
+      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
+         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
+      !#LB: 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))
+
+      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
+      t_star = dt_zu*ztmp0
+      q_star = dq_zu*ztmp0
+
+      ! 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)
+         zeta_t = zt*zeta_u/zu
+         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
+         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_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 )
+      END IF
+
+      !! ITERATION BLOCK
+      DO j_itt = 1, nb_itt
+
+         !!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
+
+         !! Update wind at zu taking into acount convection-related wind gustiness:
+         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
+         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2
+         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
+         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
+         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
+
+         !! Stability parameters:
+         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_wp), zeta_t )
+         END IF
+
+         !! Adjustment the wind at 10m (not needed in the current algo form):
+         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
+         
+         !! Roughness lengthes z0, z0t (z0q = z0t) :
+         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
+         z0    = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6)
+         ztmp1 = z0*u_star/znu_a                                           ! Re_r: roughness Reynolds number
+         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!!
+
+         !! Turbulent scales at zu :
+         ztmp0   = psi_h_coare(zeta_u)
+         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LOLO: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
+
+         t_star = dt_zu*ztmp1
+         q_star = dq_zu*ztmp1
+         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
+
+         IF( .NOT. l_zt_equal_zu ) THEN
+            ! What's need to be done if zt /= zu
+            !! Re-updating temperature and humidity at zu :
+            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
+            ztmp1 = log(zt/zu) + ztmp2
+            t_zu = t_zt - t_star/vkarmn*ztmp1
+            q_zu = q_zt - q_star/vkarmn*ztmp1
+         END IF
+
+         !! SKIN related part
+         !! -----------------
+         IF( l_use_skin ) THEN
+            !! compute transfer coefficients at zu : lolo: verifier...
+            ztmp0 = u_star/U_blk
+            Ch   = ztmp0*t_star/dt_zu
+            Ce   = ztmp0*q_star/dq_zu
+            ! Non-Solar heat flux to the ocean:
+            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
+            ztmp2 = T_s*T_s
+            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
+               &       +      rad_lw - emiss_w*stefan*ztmp2*ztmp2                           ! Net longwave flux
+            !!         => "ztmp1" is the net non-solar surface heat flux !
+            !! Updating the values of the skin temperature T_s and q_s :
+            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
+            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
+
+      ! compute transfer coefficients at zu :
+      ztmp0 = u_star/U_blk
+      Cd   = ztmp0*ztmp0
+      Ch   = ztmp0*t_star/dt_zu
+      Ce   = ztmp0*q_star/dq_zu
+      !
+      ztmp1 = zu + z0
+      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
+      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
+      Cen = Chn
+      !
+      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
+      !
+   END SUBROUTINE turb_coare3p0
+
+
+   FUNCTION alfa_charn_3p0( pwnd )
+      !!-------------------------------------------------------------------
+      !! Compute the Charnock parameter as a function of the wind speed
+      !!
+      !! (Fairall et al., 2003 p.577-578)
+      !!
+      !! Wind below 10 m/s :  alfa = 0.011
+      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
+      !! Wind greater than 18 m/s :  alfa = 0.018
+      !!
+      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p0
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) :: zw, zgt10, zgt18
+      !!-------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            zw = pwnd(ji,jj)   ! wind speed
+            !
+            ! Charnock's constant, increases with the wind :
+            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_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
+               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
+               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION alfa_charn_3p0
+
+   FUNCTION psi_m_coare( pzeta )
+      !!----------------------------------------------------------------------------------
+      !! ** Purpose: compute the universal profile stability function for momentum
+      !!             COARE 3.0, Fairall et al. 2003
+      !!             pzeta : stability paramenter, z/L where z is altitude
+      !!                     measurement and L is M-O length
+      !!       Stability function for wind speed and scalars matching Kansas and free
+      !!       convection forms with weighting f convective form, follows Fairall et
+      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
+      !!       form from Beljaars and Holtslag (1991)
+      !!
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
+      !
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            zta = pzeta(ji,jj)
+            !
+            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
+            !
+            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
+               & - 2.*ATAN(zphi_m) + 0.5*rpi
+            !
+            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
+            !
+            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+            !
+            zf = zta*zta
+            zf = zf/(1. + zf)
+            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)
+               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
+               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION psi_m_coare
+
+
+   FUNCTION psi_h_coare( pzeta )
+      !!---------------------------------------------------------------------
+      !! Universal profile stability function for temperature and humidity
+      !! COARE 3.0, Fairall et al. 2003
+      !!
+      !! pzeta : stability paramenter, z/L where z is altitude measurement
+      !!         and L is M-O length
+      !!
+      !! Stability function for wind speed and scalars matching Kansas and free
+      !! convection forms with weighting f convective form, follows Fairall et
+      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
+      !! form from Beljaars and Holtslag (1991)
+      !!
+      !! Author: L. Brodeau, June 2016 / AeroBulk
+      !!         (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            zta = pzeta(ji,jj)
+            !
+            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
+            !
+            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
+            !
+            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
+            !
+            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+            !
+            zf = zta*zta
+            zf = zf/(1. + zf)
+            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 ) &
+               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
+               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION psi_h_coare
+
+   !!======================================================================
+END MODULE sbcblk_algo_coare3p0
Index: MO/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 11283)
+++ 	(revision )
@@ -1,359 +1,0 @@
-MODULE sbcblk_algo_coare3p5
-   !!======================================================================
-   !!                       ***  MODULE  sbcblk_algo_coare3p5  ***
-   !! Computes turbulent components of surface fluxes
-   !!         according to Edson et al. 2013 (COARE v3.5) /JPO
-   !!
-   !!   * bulk transfer coefficients C_D, C_E and C_H
-   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
-   !!   * the effective bulk wind speed at 10m U_blk
-   !!   => all these are used in bulk formulas in sbcblk.F90
-   !!
-   !!    Using the bulk formulation/param. of COARE v3.5, Edson et al. 2013
-   !!
-   !!
-   !!       Routine turb_coare3p5 maintained and developed in AeroBulk
-   !!                     (https://github.com/brodeau/aerobulk/)
-   !!
-   !!            Author: Laurent Brodeau, 2016, brodeau@gmail.com
-   !!
-   !!======================================================================
-   !! History :  3.6  !  2016-02  (L.Brodeau)  Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   turb_coare3p5  : computes the bulk turbulent transfer coefficients
-   !!                   adjusts t_air and q_air from zt to zu m
-   !!                   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
-#if defined key_si3 || defined key_cice
-   USE sbc_ice         ! Surface boundary condition: ice fields
-#endif
-   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
-   USE iom             ! I/O manager library
-   USE lib_mpp         ! distribued memory computing library
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE lib_fortran     ! to use key_nosignedzero
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC ::   TURB_COARE3P5   ! called by sbcblk.F90
-
-   !                                   ! COARE own values for given constants:
-   REAL(wp), PARAMETER ::   charn0_max =   0.028   ! value above which the Charnock paramter levels off for winds > 18
-   REAL(wp), PARAMETER ::   zi0        = 600.      ! scale height of the atmospheric boundary layer...1
-   REAL(wp), PARAMETER ::   Beta0      =   1.25    ! gustiness parameter
-
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu,  &
-      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,       &
-      &                      Cdn, Chn, Cen                        )
-      !!----------------------------------------------------------------------------------
-      !!                      ***  ROUTINE  turb_coare3p5  ***
-      !!
-      !! ** 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
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
-      !!
-      !! INPUT :
-      !! -------
-      !!    *  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]
-      !!    *  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]
-      !!
-      !!
-      !! OUTPUT :
-      !! --------
-      !!    *  Cd     : drag coefficient
-      !!    *  Ch     : sensible heat coefficient
-      !!    *  Ce     : evaporation coefficient
-      !!    *  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]
-      !!
-      !!----------------------------------------------------------------------------------
-      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(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(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)
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
-      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) ::   Cdn, Chn, Cen ! neutral transfer coefficients
-      !
-      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
-         &  z0, z0t
-      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
-      !!----------------------------------------------------------------------------------
-      !
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
-
-      IF( .NOT. l_zt_equal_zu )   ALLOCATE( zeta_t(jpi,jpj) )
-
-      !! 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)  !               "
-
-      !! 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)
-
-      ztmp2   = 10000.     ! 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)
-
-      !! COARE 3.5 first guess of UN10 is U_zu
-      ztmp2 = MIN( 0.0017*U_zu - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
-      ztmp2 = MAX( ztmp2 , 0. )                       ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
-      z0     = ztmp2*u_star*u_star/grav + 0.11*znu_a/u_star
-      z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t !
-
-      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 = 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))
-
-      !! 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))
-
-      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
-      t_star = dt_zu*ztmp0
-      q_star = dq_zu*ztmp0
-
-      ! What's need to be done if zt /= zu:
-      IF( .NOT. l_zt_equal_zu ) THEN
-
-         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
-         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 )
-
-      END IF
-
-      !! ITERATION BLOCK
-      DO j_itt = 1, nb_itt
-
-         !!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]
-
-         ztmp1 = u_star*u_star   ! u*^2
-
-         !! 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
-         !!   ! 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
-         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
-
-         !! COARE 3.5: Charnock parameter is computed from the neutral wind speed at 10m: Eq. 13 (Edson al. 2013)
-         ztmp2 = u_star/vkarmn*LOG(10./z0)   ! UN10 Neutral wind at 10m!
-         ztmp2 = MIN( 0.0017*ztmp2 - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
-         ztmp2 = MAX( ztmp2 , 0. )
-
-         !! 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) )   ! COARE 3.0
-         !! Chris Fairall and Jim Edsson, private communication, March 2016 / COARE 3.5 :
-         z0t   = MIN( 1.6e-4 , 5.8E-5*ztmp1**(-0.72)) ! These thermal roughness lengths give Stanton and
-         !z0q = z0t                                   ! Dalton numbers that closely approximate COARE3.0
-
-         !! Stability parameters:
-         zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u )
-         IF( .NOT. l_zt_equal_zu ) THEN
-            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),50.0), zeta_t )
-         END IF
-
-         !! Turbulent scales at zu=10m :
-         ztmp0   = psi_h_coare(zeta_u)
-         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))
-
-         IF( .NOT. l_zt_equal_zu ) THEN
-            ! What's need to be done if zt /= zu
-            !! Re-updating temperature and humidity at zu :
-            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
-            ztmp1 = log(zt/zu) + ztmp2
-            t_zu = t_zt - t_star/vkarmn*ztmp1
-            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
-
-      END DO
-      !
-      ! compute transfer coefficients at zu :
-      ztmp0 = u_star/U_blk
-      Cd   = ztmp0*ztmp0
-      Ch   = ztmp0*t_star/dt_zu
-      Ce   = ztmp0*q_star/dq_zu
-      !
-      ztmp1 = zu + z0
-      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
-      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
-      Cen = Chn
-      !
-      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
-      !
-   END SUBROUTINE turb_coare3p5
-
-   FUNCTION psi_m_coare( pzeta )
-      !!----------------------------------------------------------------------------------
-      !! ** Purpose: compute the universal profile stability function for momentum
-      !!             COARE 3.0, Fairall et al. 2003
-      !!             pzeta : stability paramenter, z/L where z is altitude
-      !!                     measurement and L is M-O length
-      !!       Stability function for wind speed and scalars matching Kansas and free
-      !!       convection forms with weighting f convective form, follows Fairall et
-      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
-      !!       form from Beljaars and Holtslag (1991)
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zta = pzeta(ji,jj)
-            !
-            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
-            !
-            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
-               & - 2.*ATAN(zphi_m) + 0.5*rpi
-            !
-            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
-            !
-            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-            !
-            zf = zta*zta
-            zf = zf/(1. + zf)
-            zc = MIN(50., 0.35*zta)
-            zstab = 0.5 + SIGN(0.5, zta)
-            !
-            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
-               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
-               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION psi_m_coare
-
-
-   FUNCTION psi_h_coare( pzeta )
-      !!---------------------------------------------------------------------
-      !! Universal profile stability function for temperature and humidity
-      !! COARE 3.0, Fairall et al. 2003
-      !!
-      !! pzeta : stability paramenter, z/L where z is altitude measurement
-      !!         and L is M-O length
-      !!
-      !! Stability function for wind speed and scalars matching Kansas and free
-      !! convection forms with weighting f convective form, follows Fairall et
-      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
-      !! form from Beljaars and Holtslag (1991)
-      !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------
-      !!
-      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
-      !
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zta = pzeta(ji,jj)
-            !
-            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
-            !
-            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
-            !
-            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
-            !
-            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-            !
-            zf = zta*zta
-            zf = zf/(1. + zf)
-            zc = MIN(50.,0.35*zta)
-            zstab = 0.5 + SIGN(0.5, zta)
-            !
-            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
-               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
-               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION psi_h_coare
-
-   !!======================================================================
-END MODULE sbcblk_algo_coare3p5
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11284)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11284)
@@ -0,0 +1,437 @@
+MODULE sbcblk_algo_coare3p6
+   !!======================================================================
+   !!                   ***  MODULE  sbcblk_algo_coare3p6  ***
+   !!
+   !!                After Fairall et al 2018 & Edson et al 2013
+   !! 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
+   !!   * the effective bulk wind speed at 10m U_blk
+   !!   => all these are used in bulk formulas in sbcblk.F90
+   !!
+   !!       Routine turb_coare3p6 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
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   turb_coare3p6  : computes the bulk turbulent transfer coefficients
+   !!                   adjusts t_air and q_air from zt to zu m
+   !!                   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 iom             ! I/O manager library
+   USE lib_mpp         ! distribued memory computing library
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+   USE sbcwave, ONLY   :  cdn_wave ! wave module
+#if defined key_si3 || defined key_cice
+   USE sbc_ice         ! Surface boundary condition: ice fields
+#endif
+   USE lib_fortran     ! to use key_nosignedzero
+
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
+   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC ::   TURB_COARE3P6   ! called by sbcblk.F90
+
+   !                                              !! COARE own values for given constants:
+   REAL(wp), PARAMETER ::   zi0     = 600._wp     ! scale height of the atmospheric boundary layer...
+   REAL(wp), PARAMETER ::   Beta0   =   1.2_wp   ! gustiness parameter
+
+   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
+
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE turb_coare3p6( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
+      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
+      &                      Cdn, Chn, Cen,                      &
+      &                      Qsw, rad_lw, slp, Tsk_b             )
+      !!----------------------------------------------------------------------
+      !!                      ***  ROUTINE  turb_coare3p6  ***
+      !!
+      !! ** 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
+      !!                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 :
+      !! -------
+      !!    *  zt   : height for temperature and spec. hum. of air            [m]
+      !!    *  zu   : height for wind speed (usually 10m)                     [m]
+      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
+      !!    *  t_zt : potential air temperature at zt                         [K]
+      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
+      !!
+      !! 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_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]
+      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
+      !!
+      !! OUTPUT :
+      !! --------
+      !!    *  Cd     : drag coefficient
+      !!    *  Ch     : sensible heat coefficient
+      !!    *  Ce     : evaporation coefficient
+      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
+      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
+      !!    *  U_blk  : bulk wind speed at zu                                 [m/s]
+      !!
+      !!
+      !! ** Author: L. Brodeau, June 2019 / 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(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(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)
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
+      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 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(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
+      !
+      INTEGER :: j_itt
+      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
+         &  z0, z0t
+      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
+      !
+      ! Cool skin:
+      LOGICAL :: l_use_skin = .FALSE.
+      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
+      !!----------------------------------------------------------------------------------
+      !
+      ! 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
+
+      l_zt_equal_zu = .FALSE.
+      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+
+      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
+
+      !! Initialization for cool skin:
+      zsst   = T_s    ! save the bulk SST
+      IF( l_use_skin ) THEN
+         ! First guess for skin temperature:
+         IF( PRESENT(Tsk_b) ) THEN
+            T_s = Tsk_b
+         ELSE
+            T_s = T_s - 0.25     ! sst - 0.25
+         END IF
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
+      END IF
+
+      !! First guess of temperature and humidity at height zu:
+      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
+      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 - 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     = alfa_charn_3p6(U_zu)*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_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 )
+      ztmp0 = ztmp0*ztmp2
+      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
+         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
+      !#LB: 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))
+
+      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
+      t_star = dt_zu*ztmp0
+      q_star = dq_zu*ztmp0
+
+      ! 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)
+         zeta_t = zt*zeta_u/zu
+         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
+         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_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 )
+      END IF
+
+      !! ITERATION BLOCK
+      DO j_itt = 1, nb_itt
+
+         !!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
+
+         !! Update wind at zu taking into acount convection-related wind gustiness:
+         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
+         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2
+         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
+         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
+         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
+
+         !! Stability parameters:
+         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_wp), zeta_t )
+         END IF
+
+         !! Adjustment the wind at 10m (not needed in the current algo form):
+         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
+         
+         !! Roughness lengthes z0, z0t (z0q = z0t) :
+         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
+         z0    = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6)
+         ztmp1 = z0*u_star/znu_a                                           ! Re_r: roughness Reynolds number
+         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1**(-0.72_wp))   ! COARE 3.6
+
+         !! Turbulent scales at zu :
+         ztmp0   = psi_h_coare(zeta_u)
+         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LOLO: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
+
+         t_star = dt_zu*ztmp1
+         q_star = dq_zu*ztmp1
+         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
+
+         IF( .NOT. l_zt_equal_zu ) THEN
+            ! What's need to be done if zt /= zu
+            !! Re-updating temperature and humidity at zu :
+            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
+            ztmp1 = log(zt/zu) + ztmp2
+            t_zu = t_zt - t_star/vkarmn*ztmp1
+            q_zu = q_zt - q_star/vkarmn*ztmp1
+         END IF
+
+         !! SKIN related part
+         !! -----------------
+         IF( l_use_skin ) THEN
+            !! compute transfer coefficients at zu : lolo: verifier...
+            ztmp0 = u_star/U_blk
+            Ch   = ztmp0*t_star/dt_zu
+            Ce   = ztmp0*q_star/dq_zu
+            ! Non-Solar heat flux to the ocean:
+            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
+            ztmp2 = T_s*T_s
+            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
+               &       +      rad_lw - emiss_w*stefan*ztmp2*ztmp2                           ! Net longwave flux
+            !!         => "ztmp1" is the net non-solar surface heat flux !
+            !! Updating the values of the skin temperature T_s and q_s :
+            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
+            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
+
+      ! compute transfer coefficients at zu :
+      ztmp0 = u_star/U_blk
+      Cd   = ztmp0*ztmp0
+      Ch   = ztmp0*t_star/dt_zu
+      Ce   = ztmp0*q_star/dq_zu
+      !
+      ztmp1 = zu + z0
+      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
+      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
+      Cen = Chn
+      !
+      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
+      !
+   END SUBROUTINE turb_coare3p6
+
+
+   FUNCTION alfa_charn_3p6( pwnd )
+      !!-------------------------------------------------------------------
+      !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m
+      !!
+      !!  (Eq. 13 in Edson et al., 2013)
+      !!
+      !! Author: L. Brodeau, July 2019 / AeroBulk  (https://github.com/brodeau/aerobulk/)
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! neutral wind speed at 10m
+      !
+      REAL(wp), PARAMETER :: charn0_max = 0.028  !: value above which the Charnock parameter levels off for winds > 18 m/s
+      !!-------------------------------------------------------------------
+      !
+      alfa_charn_3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp )
+      !
+   END FUNCTION alfa_charn_3p6
+
+   FUNCTION psi_m_coare( pzeta )
+      !!----------------------------------------------------------------------------------
+      !! ** Purpose: compute the universal profile stability function for momentum
+      !!             COARE 3.0, Fairall et al. 2003
+      !!             pzeta : stability paramenter, z/L where z is altitude
+      !!                     measurement and L is M-O length
+      !!       Stability function for wind speed and scalars matching Kansas and free
+      !!       convection forms with weighting f convective form, follows Fairall et
+      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
+      !!       form from Beljaars and Holtslag (1991)
+      !!
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
+      !
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            zta = pzeta(ji,jj)
+            !
+            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
+            !
+            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
+               & - 2.*ATAN(zphi_m) + 0.5*rpi
+            !
+            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
+            !
+            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+            !
+            zf = zta*zta
+            zf = zf/(1. + zf)
+            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)
+               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
+               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION psi_m_coare
+
+
+   FUNCTION psi_h_coare( pzeta )
+      !!---------------------------------------------------------------------
+      !! Universal profile stability function for temperature and humidity
+      !! COARE 3.0, Fairall et al. 2003
+      !!
+      !! pzeta : stability paramenter, z/L where z is altitude measurement
+      !!         and L is M-O length
+      !!
+      !! Stability function for wind speed and scalars matching Kansas and free
+      !! convection forms with weighting f convective form, follows Fairall et
+      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
+      !! form from Beljaars and Holtslag (1991)
+      !!
+      !! Author: L. Brodeau, June 2016 / AeroBulk
+      !!         (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            zta = pzeta(ji,jj)
+            !
+            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
+            !
+            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
+            !
+            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
+            !
+            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+            !
+            zf = zta*zta
+            zf = zf/(1. + zf)
+            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 ) &
+               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
+               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION psi_h_coare
+
+   !!======================================================================
+END MODULE sbcblk_algo_coare3p6
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11283)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11284)
@@ -17,5 +17,6 @@
    !!   Ri_bulk       : bulk Richardson number aka BRN
    !!   q_sat         : saturation humidity as a function of SLP and temperature
-
+   !!   q_air_rh      : specific humidity as a function of RH, t_air and SLP
+   
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -37,4 +38,5 @@
    PUBLIC Ri_bulk
    PUBLIC q_sat
+   PUBLIC q_air_rh
 
    !!----------------------------------------------------------------------
@@ -115,5 +117,5 @@
       !!                           ***  FUNCTION L_vap  ***
       !!
-      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
+      !! ** Purpose : Compute the latent heat of vaporization of water out of temperature
       !!
       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
@@ -273,4 +275,33 @@
 
 
+   FUNCTION e_sat_sclr( ptak, pslp )
+      !!----------------------------------------------------------------------------------
+      !!                   ***  FUNCTION e_sat_sclr  ***
+      !!                  < SCALAR argument version >
+      !! ** Purpose : water vapor at saturation in [Pa]
+      !!              Based on accurate estimate by Goff, 1957
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!
+      !!    Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here
+      !!----------------------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K]
+      REAL(wp), INTENT(in) ::   pslp    ! sea level atmospheric pressure   [Pa]
+      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg]
+      !
+      REAL(wp) ::   zta, ztmp   ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
+      ztmp = rt0 / zta
+      !
+      ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957)
+      e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        &
+         &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  &
+         &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
+      !
+   END FUNCTION e_sat_sclr
+
+
    FUNCTION q_sat( ptak, pslp )
       !!----------------------------------------------------------------------------------
@@ -288,19 +319,13 @@
       !
       INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::   zta, ze_sat, ztmp   ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zta = MAX( ptak(ji,jj) , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
-            ztmp = rt0 / zta
-            !
-            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
-            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        &
-               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/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]
+      REAL(wp) ::   ze_sat   ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            ze_sat =  e_sat_sclr( ptak(ji,jj), pslp(ji,jj) )
+            !
+            q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat )
             !
          END DO
@@ -308,4 +333,28 @@
       !
    END FUNCTION q_sat
+
+   FUNCTION q_air_rh(prha, ptak, pslp)
+      !!----------------------------------------------------------------------------------
+      !! Specific humidity of air out of Relative Humidity
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             :: q_air_rh
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha        !: relative humidity      [fraction, not %!!!]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak        !: air temperature        [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp        !: atmospheric pressure   [Pa]
+      !
+      INTEGER  ::   ji, jj      ! dummy loop indices
+      REAL(wp) ::   ze      ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj))
+            q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze)
+         END DO
+      END DO
+      !
+   END FUNCTION q_air_rh
 
 
