Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11182)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11209)
@@ -31,5 +31,5 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
 #endif
-   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
    USE in_out_manager ! I/O manager
    USE iom            ! I/O manager library
@@ -149,11 +149,12 @@
       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
 
-      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
-
-      !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2):
+      ztmp2 = Ri_bulk( zu, sst, t_zu, ssq, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
+      !ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
+
+      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
       ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
       ztmp0 = ztmp0*ztmp2
-      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  Ribu < 0
-         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  Ribu > 0
+      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  BRN < 0
+         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  BRN > 0
       !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
 
@@ -285,37 +286,4 @@
    END FUNCTION alfa_charn
 
-
-   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
-      !!------------------------------------------------------------------------
-      !!
-      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
-      !!  specific humidity, and frictional scales u*, t* and q*
-      !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://github.com/brodeau/aerobulk/)
-      !!------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
-         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
-         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::     zqa          ! local scalar
-      !!-------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zqa = (1. + rctv0*pqa(ji,jj))
-            !
-            One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
-               &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION One_on_L
-
-
    FUNCTION psi_m_coare( pzeta )
       !!----------------------------------------------------------------------------------
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90	(revision 11182)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90	(revision 11209)
@@ -35,5 +35,5 @@
    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
@@ -255,37 +255,4 @@
    END SUBROUTINE turb_coare3p5
 
-
-
-   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
-      !!------------------------------------------------------------------------
-      !!
-      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
-      !!  specific humidity, and frictional scales u*, t* and q*
-      !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://github.com/brodeau/aerobulk/)
-      !!------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
-         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
-         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::     zqa          ! local scalar
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zqa = (1. + rctv0*pqa(ji,jj))
-            !
-            One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
-               &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION One_on_L
-
-
    FUNCTION psi_m_coare( pzeta )
       !!----------------------------------------------------------------------------------
@@ -389,28 +356,4 @@
    END FUNCTION psi_h_coare
 
-
-   FUNCTION visc_air( ptak )
-      !!---------------------------------------------------------------------
-      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
-      !!
-      !! Author: L. Brodeau, june 2016 / AeroBulk
-      !!         (https://github.com/brodeau/aerobulk/)
-      !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature   [K]
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::   ztc, ztc2      ! local scalar
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
-            ztc2 = ztc*ztc
-            visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc)
-         END DO
-      END DO
-      !
-   END FUNCTION visc_air
-
    !!======================================================================
 END MODULE sbcblk_algo_coare3p5
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11182)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11209)
@@ -38,5 +38,5 @@
 
    USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbcblk_phy   !LB: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
    USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF)
 
@@ -72,8 +72,4 @@
    !                            !:    => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions)
    !                            !:  => Otherwize use "nb_itt_wl = 10"
-
-
-
-
    !!----------------------------------------------------------------------
 CONTAINS
@@ -213,11 +209,11 @@
       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
 
-      ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk )   ! Ribu = Bulk Richardson number
-
-      !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2):
+      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
+
+      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
       ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
       func_m = ztmp0*ztmp2 ! temporary array !!
-      func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  Ribu < 0 ! temporary array !!! func_h == zeta_u
-         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  Ribu > 0
+      func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u
+         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))           !  BRN > 0
       !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
 
@@ -247,7 +243,5 @@
 
       !! First guess of inverse of Monin-Obukov length (1/L) :
-      ztmp0 = (1. + rctv0*q_zu)  ! the factor to apply to temp. to get virt. temp...
-      Linv  =  grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / MAX( u_star*u_star * t_zu*ztmp0 , 1.E-9 ) ! #LOLO
-      Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
+      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
 
       !! Functions such as  u* = U_blk*vkarmn/func_m
@@ -260,5 +254,5 @@
 
          !! Bulk Richardson Number at z=zu (Eq. 3.25)
-         ztmp0 = Ri_bulk(zu, t_zu, dt_zu, q_zu, dq_zu, U_blk)
+         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
 
          !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :
@@ -278,5 +272,5 @@
          z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
          z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
-         
+
          !! Update wind at 10m taking into acount convection-related wind gustiness:
          !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8
@@ -447,25 +441,4 @@
 
 
-   FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub )
-      !!----------------------------------------------------------------------------------
-      !! Bulk Richardson number (Eq. 3.25 IFS doc)
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) ::   Ri_bulk   !
-      !
-      REAL(wp)                    , INTENT(in) ::   pz    ! height above the sea        [m]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptz   ! air temperature at pz m     [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdt   ! ptz - sst                   [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqz   ! air temperature at pz m [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdq   ! pqz - ssq               [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pub   ! bulk wind speed           [m/s]
-      !!----------------------------------------------------------------------------------
-      !
-      Ri_bulk =   grav*pz/(pub*pub)   &
-         &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(rCp_dry + rCp_vap*pqz))) &
-         &          + rctv0*pdq )
-      !
-   END FUNCTION Ri_bulk
 
 
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11182)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11209)
@@ -38,4 +38,5 @@
    USE lib_fortran     ! to use key_nosignedzero
 
+   USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
 
    IMPLICIT NONE
@@ -93,5 +94,5 @@
       !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
       !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
-      !!    *  U_blk  : bulk wind at 10m                                      [m/s]
+      !!    *  U_blk  : bulk wind speed at 10m                                [m/s]
       !!----------------------------------------------------------------------------------
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
@@ -125,9 +126,9 @@
       IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
 
-      U_blk = MAX( 0.5 , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
+      U_blk = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
 
       !! First guess of stability:
-      ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt
-      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
+      ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt
+      stab  = 0.5_wp + sign(0.5_wp,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
 
       !! Neutral coefficients at 10m:
@@ -159,15 +160,14 @@
 
          ! Updating turbulent scales :   (L&Y 2004 eq. (7))
-         ztmp1  = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd))
-         ztmp2  = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd))
-
-         ztmp0 = 1. + rctv0*q_zu      ! multiply this with t and you have the virtual temperature
+         ztmp0 = stab*U_blk       ! u*       (stab == SQRT(Cd))
+         ztmp1 = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd))
+         ztmp2 = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd))
 
          ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
-         ztmp0 =  (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk)
-         !                                                      ( Cd*U_blk*U_blk is U*^2 at zu )
-
+         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
+         
          !! Stability parameters :
-         zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )
+         zeta_u   = zu*ztmp0
+         zeta_u = sign( min(abs(zeta_u),10.0_wp), zeta_u )
          zpsi_h_u = psi_h( zeta_u )
 
@@ -175,9 +175,10 @@
          IF( .NOT. l_zt_equal_zu ) THEN
             !! Array 'stab' is free for the moment so using it to store 'zeta_t'
-            stab = zt*ztmp0 ;  stab = SIGN( MIN(ABS(stab),10.0), stab )  ! Temporaty array stab == zeta_t !!!
+            stab = zt*ztmp0
+            stab = SIGN( MIN(ABS(stab),10.0_wp), stab )  ! Temporaty array stab == zeta_t !!!
             stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again!
             t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b)
             q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c)
-            q_zu = max(0., q_zu)
+            q_zu = max(0._wp, q_zu)
          END IF
 
@@ -194,32 +195,32 @@
 
          ELSE
-            ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
-            !   In very rare low-wind conditions, the old way of estimating the
-            !   neutral wind speed at 10m leads to a negative value that causes the code
-            !   to crash. To prevent this a threshold of 0.25m/s is imposed.
-            ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))
-            ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
-            Cdn(:,:) = ztmp0
-            sqrt_Cd_n10 = sqrt(ztmp0)
-
-            stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability
-            Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10)
-            Chn(:,:) = Cx_n10
-
-            !! Update of transfer coefficients:
-            ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
-            Cd      = ztmp0 / ( ztmp1*ztmp1 )
-            stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
-
-            ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
-            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
-            ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
-            Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
-
-            Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
-            Cen(:,:) = Cx_n10
-            ztmp1 = 1. + Cx_n10*ztmp0
-            Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
-            ENDIF
+         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
+         !   In very rare low-wind conditions, the old way of estimating the
+         !   neutral wind speed at 10m leads to a negative value that causes the code
+         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
+         ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))
+         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
+         Cdn(:,:) = ztmp0
+         sqrt_Cd_n10 = sqrt(ztmp0)
+
+         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability
+         Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10)
+         Chn(:,:) = Cx_n10
+
+         !! Update of transfer coefficients:
+         ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
+         Cd      = ztmp0 / ( ztmp1*ztmp1 )
+         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
+
+         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
+         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
+         ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
+         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
+
+         Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
+         Cen(:,:) = Cx_n10
+         ztmp1 = 1. + Cx_n10*ztmp0
+         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
+         ENDIF
          !
       END DO
@@ -252,5 +253,5 @@
             !
             ! When wind speed > 33 m/s => Cyclone conditions => special treatment
-            zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) )   ! If pw10 < 33. => 0, else => 1
+            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
             !
             cd_neutral_10m(ji,jj) = 1.e-3 * ( &
@@ -258,5 +259,5 @@
                &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s
             !
-            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6)
+            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp)
             !
          END DO
@@ -285,12 +286,12 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )
-            zx2 = MAX ( zx2 , 1. )
+            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
+            zx2 = MAX( zx2 , 1._wp )
             zx  = SQRT( zx2 )
-            zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) )
-            !
-            psi_m(ji,jj) =        zstab  * (-5.*pzeta(ji,jj))       &          ! Stable
-               &          + (1. - zstab) * (2.*LOG((1. + zx)*0.5)   &          ! Unstable
-               &               + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5)  !    "
+            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
+            !
+            psi_m(ji,jj) =        zstab  * (-5._wp*pzeta(ji,jj))       &          ! Stable
+               &          + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp)   &          ! Unstable
+               &               + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp)  !    "
             !
          END DO
@@ -319,10 +320,10 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )
-            zx2 = MAX ( zx2 , 1. )
-            zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) )
-            !
-            psi_h(ji,jj) =         zstab  * (-5.*pzeta(ji,jj))        &  ! Stable
-               &           + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5 ))   ! Unstable
+            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
+            zx2 = MAX( zx2 , 1._wp )
+            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
+            !
+            psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable
+               &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable
             !
          END DO
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 11182)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11209)
@@ -8,11 +8,14 @@
    !!----------------------------------------------------------------------
 
+   !!   virt_temp     : virtual (aka sensible) temperature (potential or absolute)
    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
+   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature
+   !!   L_vap         : latent heat of vaporization of water as a function of temperature
    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
+   !!   gamma_moist   : adiabatic lapse-rate of moist air
+   !!   One_on_L      : 1. / ( Monin-Obukhov length )
+   !!   Ri_bulk       : bulk Richardson number aka BRN
    !!   q_sat         : saturation humidity as a function of SLP and temperature
-   !!   gamma_moist   : 
-   !!   L_vap         : latent heat of vaporization of water as a function of temperature
-   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature
-   
+
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -21,10 +24,17 @@
    PRIVATE
 
+   INTERFACE gamma_moist
+      MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr
+   END INTERFACE gamma_moist
+
+   PUBLIC virt_temp
    PUBLIC rho_air
+   PUBLIC visc_air
+   PUBLIC L_vap
    PUBLIC cp_air
+   PUBLIC gamma_moist
+   PUBLIC One_on_L
+   PUBLIC Ri_bulk
    PUBLIC q_sat
-   PUBLIC gamma_moist
-   PUBLIC L_vap
-   PUBLIC visc_air
 
    !!----------------------------------------------------------------------
@@ -35,4 +45,29 @@
 CONTAINS
 
+   FUNCTION virt_temp( pta, pqa )
+      !!------------------------------------------------------------------------
+      !!
+      !! Compute the (absolute/potential) virtual temperature, knowing the
+      !! (absolute/potential) temperature and specific humidity
+      !!
+      !! If input temperature is absolute then output vitual temperature is absolute
+      !! If input temperature is potential then output vitual temperature is potential
+      !!
+      !! Author: L. Brodeau, June 2019 / AeroBulk
+      !!         (https://github.com/brodeau/aerobulk/)
+      !!------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp         !: 1./(Monin Obukhov length) [m^-1]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta,  &  !: absolute or potetntial air temperature [K]
+         &                                        pqa      !: specific humidity of air   [kg/kg]
+      !!-------------------------------------------------------------------
+      !
+      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
+      !!
+      !! This is exactly the same sing that:
+      !! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa))
+      !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
+      !
+   END FUNCTION virt_temp
+
    FUNCTION rho_air( ptak, pqa, pslp )
       !!-------------------------------------------------------------------------------
@@ -41,5 +76,5 @@
       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
       !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!-------------------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
@@ -53,4 +88,42 @@
    END FUNCTION rho_air
 
+   FUNCTION visc_air(ptak)
+      !!----------------------------------------------------------------------------------
+      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
+      !
+      INTEGER  ::   ji, jj      ! dummy loop indices
+      REAL(wp) ::   ztc, ztc2   ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
+            ztc2 = ztc*ztc
+            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
+         END DO
+      END DO
+      !
+   END FUNCTION visc_air
+
+   FUNCTION L_vap( psst )
+      !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION L_vap  ***
+      !!
+      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
+      !!----------------------------------------------------------------------------------
+      !
+      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
+      !
+   END FUNCTION L_vap
 
    FUNCTION cp_air( pqa )
@@ -69,4 +142,128 @@
       !
    END FUNCTION cp_air
+
+   FUNCTION gamma_moist_vctr( ptak, pqa )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION gamma_moist_vctr  ***
+      !!
+      !! ** Purpose : Compute the moist adiabatic lapse-rate.
+      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
+      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist_vctr   ! moist adiabatic lapse-rate
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zta = MAX( ptak(ji,jj),  180._wp) ! prevents screw-up over masked regions where field == 0.
+            zqa = MAX( pqa(ji,jj),  1.E-6_wp) !    "                   "                     "
+            !
+            zwa = zqa / (1. - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
+            ziRT = 1._wp/(R_dry*zta)    ! 1/RT
+            gamma_moist_vctr(ji,jj) = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
+         END DO
+      END DO
+      !
+   END FUNCTION gamma_moist_vctr
+
+   FUNCTION gamma_moist_sclr( ptak, pqa )
+      !!----------------------------------------------------------------------------------
+      !! ** Purpose : Compute the moist adiabatic lapse-rate.
+      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
+      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
+      !!
+      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp)             :: gamma_moist_sclr
+      REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg)
+      !
+      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar
+      !!----------------------------------------------------------------------------------
+      zta = MAX( ptak,  180._wp) ! prevents screw-up over masked regions where field == 0.
+      zqa = MAX( pqa,  1.E-6_wp) !    "                   "                     "
+      !!
+      zwa = zqa / (1._wp - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
+      ziRT = 1._wp / (R_dry*zta)    ! 1/RT
+      gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta )
+      !!
+   END FUNCTION gamma_moist_sclr
+
+   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
+      !!------------------------------------------------------------------------
+      !!
+      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
+      !!  specific humidity, and frictional scales u*, t* and q*
+      !!
+      !! Author: L. Brodeau, June 2016 / AeroBulk
+      !!         (https://github.com/brodeau/aerobulk/)
+      !!------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
+         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
+         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) ::     zqa          ! local scalar
+      !!-------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            zqa = (1._wp + rctv0*pqa(ji,jj))
+            !
+            One_on_L(ji,jj) = grav*vkarmn*(pts(ji,jj) + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
+               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
+            !
+         END DO
+      END DO
+      !
+      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
+      !
+   END FUNCTION One_on_L
+
+   FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub )
+      !!----------------------------------------------------------------------------------
+      !! Bulk Richardson number according to "wide-spread equation"...
+      !!
+      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk
+      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst  ! SST                                   [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub   ! bulk wind speed                     [m/s]
+      !
+      INTEGER  ::   ji, jj                                ! dummy loop indices
+      REAL(wp) ::   zqa, zta, zgamma, zdth_v, ztv, zsstv  ! local scalars
+      !!-------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer...
+            zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(ptha(ji,jj),zqa)*pz ) ! ~ mean absolute temperature of air within the layer
+            zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer
+            zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer
+            !
+            zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!)
+            !
+            zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature"
+            !
+            ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer
+            !
+            Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk
+            !
+         END DO
+      END DO
+   END FUNCTION Ri_bulk
 
 
@@ -108,75 +305,4 @@
 
 
-   FUNCTION gamma_moist( ptak, pqa )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION gamma_moist  ***
-      !!
-      !! ** Purpose : Compute the moist adiabatic lapse-rate.
-      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
-      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) :: zrv, ziRT        ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
-            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
-            gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( rCp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )
-         END DO
-      END DO
-      !
-   END FUNCTION gamma_moist
-
-
-   FUNCTION L_vap( psst )
-      !!---------------------------------------------------------------------------------
-      !!                           ***  FUNCTION L_vap  ***
-      !!
-      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
-      !!----------------------------------------------------------------------------------
-      !
-      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
-      !
-   END FUNCTION L_vap
-
-
-
-   FUNCTION visc_air(ptak)
-      !!----------------------------------------------------------------------------------
-      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   !
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
-      !
-      INTEGER  ::   ji, jj      ! dummy loop indices
-      REAL(wp) ::   ztc, ztc2   ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
-            ztc2 = ztc*ztc
-            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
-         END DO
-      END DO
-      !
-   END FUNCTION visc_air
-   
-
    !!======================================================================
 END MODULE sbcblk_phy
