Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90	(revision 11181)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90	(revision 11182)
@@ -114,4 +114,5 @@
    !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rhoa              !: air density at "rn_zu" m above the sea       [kg/m3] !LB
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr               !: sea heat flux:     solar                     [W/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2]
@@ -173,5 +174,5 @@
       !
       ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) ,     &
-         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , STAT=ierr(1) ) 
+         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) ) 
          !
       ALLOCATE( qns_tot(jpi,jpj) , qns  (jpi,jpj) , qns_b(jpi,jpj),        &
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 11181)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11182)
@@ -56,5 +56,5 @@
    USE prtctl         ! Print control
 
-   USE sbcblk_phymbl  !LB: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_phy     !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc...
 
    
@@ -109,5 +109,5 @@
    !LB:
    LOGICAL  ::   ln_skin        ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
-   !LB: => sbc_oce.F90: REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk                       ! sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius]
+   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature if .true. (specific humidity espected if .false.) !LB
    !LB.
    
@@ -175,5 +175,5 @@
          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
-         &                 ln_skin                                                        ! cool-skin / warm-layer param. !LB
+         &                 ln_skin, ln_humi_dpt                                                        ! cool-skin / warm-layer param. !LB
       !!---------------------------------------------------------------------
       !
@@ -292,4 +292,6 @@
          WRITE(numout,*)
          WRITE(numout,*) '      use cool-skin/warm-layer parameterization (SSST)     ln_skin     = ', ln_skin !LB
+         WRITE(numout,*)
+         WRITE(numout,*) '      air humidity read in files is dew-point temperature? ln_humi_dpt = ', ln_humi_dpt !LB
          !
       ENDIF
@@ -346,5 +348,11 @@
          ENDIF 
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
-         qatm_ice(:,:)    = sf(jp_humi)%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
+         !LB.
          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
          sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
@@ -386,5 +394,5 @@
       REAL(wp) ::   zztmp                ! local variable
       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst  [kg/kg]
       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes
       REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation
@@ -392,5 +400,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
+      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg] !LB
       !!---------------------------------------------------------------------
       !
@@ -443,17 +451,28 @@
       ! ... specific humidity at SST and IST tmask(
       zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
-      !!
+
+      !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)
+         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
+      !LB.
+
       !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
       !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
       !!    (since reanalysis products provide T at z, not theta !)
-      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
+      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt
 
       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
       !
-      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
+      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
+      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( 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, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
+      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
 
@@ -463,7 +482,7 @@
             !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL optional arrays!!!'
             !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => BEFORE ZST(40:50,30) =', zst(40:50,30)
-            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF         
-               &                                           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))
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF         
+               &                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))
             !LB: "zst" and "zsq" have been updated with skin temp. !!! (from bulk SST)...
             zst(:,:) = zst(:,:)*tmask(:,:,1)
@@ -472,5 +491,5 @@
          ELSE
             IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
-            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF         
+            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF         
                &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
          END IF
@@ -487,8 +506,8 @@
 
       !                          ! Compute true air density :
-      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
-         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
+      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove rhoa*grav*rn_zu from SLP...)
+         rhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
       ELSE                                      ! At zt:
-         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+         rhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), zqair(:,:), sf(jp_slp)%fnow(:,:,1) )
       END IF
 
@@ -498,5 +517,5 @@
       DO jj = 1, jpj             ! tau module, i and j component
          DO ji = 1, jpi
-            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
+            zztmp = rhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
             taum  (ji,jj) = zztmp * wndm  (ji,jj)
             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
@@ -527,15 +546,15 @@
 
       ! zqla used as temporary array, for rho*U (common term of bulk formulae):
-      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
+      zqla(:,:) = rhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
 
       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
          !! q_air and t_air are given at 10m (wind reference height)
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
+         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - zqair(:,:)) ) ! Evaporation, using bulk wind speed
+         zqsb (:,:) = cp_air(zqair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
       ELSE
          !! q_air and t_air are not given at 10m (wind reference height)
          ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
+         zqsb (:,:) = cp_air(zqair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
       ENDIF
 
@@ -584,5 +603,5 @@
       !
       !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN
-         CALL iom_put( "rho_air" ,   zrhoa )                ! output air density (kg/m^3) !#LB
+         CALL iom_put( "rho_air" ,   rhoa )                 ! output air density (kg/m^3) !#LB
          CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
          CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
@@ -638,5 +657,4 @@
       REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
       REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau)
       !!---------------------------------------------------------------------
       !
@@ -673,6 +691,4 @@
 
       ! local scalars ( place there for vector optimisation purposes)
-      ! Computing density of air! Way denser that 1.2 over sea-ice !!!
-      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
 
       !!gm brutal....
@@ -687,7 +703,7 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
-            utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
+            utau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
                &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
-            vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
+            vtau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
                &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
          END DO
@@ -732,5 +748,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
       !!---------------------------------------------------------------------
       !
@@ -738,6 +754,16 @@
       zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
       !
-      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
-      !
+
+      !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)
+         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
+      !LB.
+
       zztmp = 1. / ( 1. - albo )
       WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
@@ -766,8 +792,8 @@
                ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
                ! Sensible Heat
-               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
+               z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
                ! Latent Heat
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
-                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
+                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / rhoa(ji,jj) - zqair(ji,jj) ) )
                ! Latent heat sensitivity for ice (Dqla/Dt)
                IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
@@ -779,5 +805,5 @@
 
                ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
-               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
+               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
 
                ! ----------------------------!
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 11181)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90	(revision 11182)
@@ -31,5 +31,5 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
 #endif
-   USE sbcblk_phymbl  !LB: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc...
    USE in_out_manager ! I/O manager
    USE iom            ! I/O manager library
@@ -104,5 +104,5 @@
       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, &
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 11181)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11182)
@@ -38,5 +38,5 @@
 
    USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbcblk_phymbl   !LB: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_phy   !LB: all thermodynamics functions, rho_air, q_sat, etc...
    USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF)
 
@@ -265,5 +265,5 @@
          Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
          !! Note: it is slightly different that the L we would get with the usual
-         !! expression, as in coare algorithm or in 'mod_phymbl.f90' (One_on_L_MO())
+         !! expression, as in coare algorithm or in 'mod_phy.f90' (One_on_L_MO())
          Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
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 11182)
@@ -0,0 +1,182 @@
+MODULE sbcblk_phy
+   !!======================================================================
+   !!                       ***  MODULE  sbcblk_phy  ***
+   !! A set of functions to compute air themodynamics parameters
+   !!                     needed by Aerodynamic Bulk Formulas
+   !!=====================================================================
+   !!            4.0  !  2019 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
+   !!----------------------------------------------------------------------
+
+   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
+   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
+   !!   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
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC rho_air
+   PUBLIC cp_air
+   PUBLIC q_sat
+   PUBLIC gamma_moist
+   PUBLIC L_vap
+   PUBLIC visc_air
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   FUNCTION rho_air( ptak, pqa, pslp )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION rho_air  ***
+      !!
+      !! ** 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/)
+      !!-------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
+      !!-------------------------------------------------------------------------------
+      !
+      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
+      !
+   END FUNCTION rho_air
+
+
+   FUNCTION cp_air( pqa )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air  ***
+      !!
+      !! ** Purpose : Compute specific heat (Cp) of moist air
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
+      !!-------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
+      !!-------------------------------------------------------------------------------
+      !
+      cp_air = rCp_dry + rCp_vap * pqa
+      !
+   END FUNCTION cp_air
+
+
+   FUNCTION q_sat( ptak, pslp )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION q_sat  ***
+      !!
+      !! ** Purpose : Specific humidity at saturation in [kg/kg]
+      !!              Based on accurate estimate of "e_sat"
+      !!              aka saturation water vapor (Goff, 1957)
+      !!
+      !! ** 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) ::   pslp    ! sea level atmospheric pressure       [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
+      !
+      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]
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION q_sat
+
+
+   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
Index: MO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phymbl.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phymbl.F90	(revision 11181)
+++ 	(revision )
@@ -1,193 +1,0 @@
-MODULE sbcblk_phymbl
-   !!======================================================================
-   !!                       ***  MODULE  sbcblk_phymbl  ***
-   !! A set of functions to compute air themodynamics parameters
-   !!                     needed by Aerodynamic Bulk Formulas
-   !!=====================================================================
-   !!            4.0  !  2019 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
-   !!----------------------------------------------------------------------
-
-   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
-   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
-   !!   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 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 sbcdcy         ! surface boundary condition: diurnal cycle
-   !USE sbcwave , ONLY :   cdn_wave ! wave module
-   !USE sbc_ice        ! Surface boundary condition: ice fields
-   USE lib_fortran    ! to use key_nosignedzero
-   !
-   !USE iom            ! I/O manager library
-   !USE in_out_manager ! I/O manager
-   USE lib_mpp        ! distribued memory computing library
-   !USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   !USE prtctl         ! Print control
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC rho_air
-   PUBLIC cp_air
-   PUBLIC q_sat
-   PUBLIC gamma_moist
-   PUBLIC L_vap
-   PUBLIC visc_air
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   FUNCTION rho_air( ptak, pqa, pslp )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION rho_air  ***
-      !!
-      !! ** 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/)
-      !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
-      !!-------------------------------------------------------------------------------
-      !
-      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
-      !
-   END FUNCTION rho_air
-
-
-   FUNCTION cp_air( pqa )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION cp_air  ***
-      !!
-      !! ** Purpose : Compute specific heat (Cp) of moist air
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
-      !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
-      !!-------------------------------------------------------------------------------
-      !
-      cp_air = rCp_dry + rCp_vap * pqa
-      !
-   END FUNCTION cp_air
-
-
-   FUNCTION q_sat( ptak, pslp )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION q_sat  ***
-      !!
-      !! ** Purpose : Specific humidity at saturation in [kg/kg]
-      !!              Based on accurate estimate of "e_sat"
-      !!              aka saturation water vapor (Goff, 1957)
-      !!
-      !! ** 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) ::   pslp    ! sea level atmospheric pressure       [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::   ze_sat, ztmp   ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            ztmp = rt0 / ptak(ji,jj)
-            !
-            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
-            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
-               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/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]
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION q_sat
-
-
-   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_phymbl
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin.F90	(revision 11181)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin.F90	(revision 11182)
@@ -25,5 +25,5 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
 
-   USE sbcblk_phymbl  !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
 
    USE lib_mpp         ! distribued memory computing library
