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