module sens_heat_rain_m

  implicit none

contains

  elemental real function sens_heat_rain(rain, t, q, rhoa, xlv, t_int, p)

    ! Computes heat flux due to rainfall, in W m-2, positive upward.

    use const, only: cpa, cpw, rgas
#ifndef IN_LMDZ
    use const, only: eps_w
#endif
    use esat_m, only: esat

#ifdef IN_LMDZ
    include "YOMCST.h"
    ! for eps_w
#endif

    real, intent(in):: rain ! rain mass flux, in kg m-2 s-1
    real, intent(in):: t ! air temperature at 10 m, in K
    real, intent(in):: q ! specific humidity at 10 m
    real, intent(in):: rhoa ! density of moist air  (kg / m3)
    real, intent(in):: xlv ! latent heat of evaporation (J / kg)
    real, intent(in):: t_int ! interface temperature, in K
    real, intent(in):: p ! surface pressure, in Pa

    ! Local:
    
    real es ! saturation pressure of wator vapor, in Pa
    real alfac ! wet bulb factor
    real dwat ! water vapour diffusivity
    real dtmp ! heat diffusivity
    real q_int ! specific (saturation) humidity at ocean interface
    real t_celsius ! air temperature at 10 m, in Celsius degrees

    real wetc
    ! derivative of saturated mass fraction of water vapor with
    ! respect to temperature, at constant total pressure

    !---------------------------------------------------------------------

    es = esat(t_int, p) * 0.98 ! reduced for salinity, Kraus 1972 page 46
    q_int = eps_w * (es / (p - (1. - eps_w) * es))
    wetc = eps_w * xlv * q_int / (rgas * t_int**2)
    dwat = 2.11e-5 * (t / 273.15)**1.94
    t_celsius = t - 273.15
    dtmp = (1. + 3.309e-3 * t_celsius - 1.44e-6 * t_celsius**2) * 0.02411 &
         / (rhoa * cpa)

    ! Gosnell 1995 k0991, equation (11):
    alfac =  1. / (1. + (wetc * xlv * dwat) / (cpa * dtmp))

    ! Gosnell 1995 k0991, equation (12):
    sens_heat_rain =  rain * alfac * cpw * (t_int - t + (q_int - q) * xlv / cpa)

  end function sens_heat_rain
  
end module sens_heat_rain_m