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 11266)
+++ 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
 
 
