module Microlayer_m

  Implicit none

contains

  subroutine Microlayer(dter, dser, tkt, tks, hlb, tau, s_subskin, al, &
       xlv, taur, rf, rain, qcol)

    ! H. Bellenger 2016

    use const, only: beta, cpw, grav, rhow
    use fv_m, only: fv

    real, intent(out):: dter(:)
    ! Temperature variation in the diffusive microlayer, that is
    ! ocean-air interface temperature minus subskin temperature. In K.

    real, intent(out):: dser(:)
    ! Salinity variation in the diffusive microlayer, that is ocean-air
    ! interface salinity minus subskin salinity. In ppt.

    real, intent(inout):: tkt(:)
    ! thickness of cool skin (microlayer), in m

    real, intent(inout):: tks(:)
    ! thickness of mass diffusion layer (microlayer), in m

    real, intent(in):: hlb(:)
    ! latent heat flux at the surface, positive upward (W m-2)

    real, intent(in):: tau(:) ! wind stress, turbulent part only, in Pa
    real, intent(in):: s_subskin(:) ! subskin salinity, in ppt
    real, intent(in):: al(:) ! water thermal expansion coefficient (in K-1)
    real, intent(in):: xlv(:) ! latent heat of evaporation (J/kg)
    real, intent(in):: taur(:) ! momentum flux due to rainfall, in Pa

    real, intent(in):: rf(:)
    ! sensible heat flux at the surface due to rainfall, in W m-2

    real, intent(in):: rain(:) ! rain mass flux, in kg m-2 s-1

    real, intent(in):: qcol(:)
    ! net flux at the surface, without sensible heat flux due to rain, in W m-2

    ! Local:

    real, dimension(size(qcol)):: usrk, usrct, usrcs, alq
    real xlamx(size(qcol)) ! Saunders coefficient
    real, parameter:: visw = 1e-6
    real, parameter:: tcw = 0.6 ! thermal conductivity of water

    real, parameter:: mu = 0.0129e-7 ! in m2 / s
    ! molecular salinity diffusivity, Kraus and Businger, page 47

    real, parameter:: kappa = 1.49e-7 ! thermal diffusivity, in m2 / s

    real, parameter:: afk = 4e-4
    real, parameter:: bfk = 1.3
    ! a and b coefficient for the power function fitting the TKE flux 
    ! carried by rain:  Fk = a * R**b, derived form the exact solution
    ! of Soloviev and Lukas 2006 (Schlussel et al 1997, Craeye and
    ! Schlussel 1998)

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

    alq = al * (qcol + rf * (1 - fV(tkt, rain))) - beta * s_subskin * cpw &
         * (hlb / xlv - rain * (1 - fV(tks, rain)))

    usrk = (afk / rhow)**(1. / 3.) * (rain * 3600.)**(bfk / 3.)
    ! Equivalent friction velocity due to the TKE input by the penetrating
    ! raindrops Fk

    ! Friction velocities in the air:
    usrct = sqrt((tau + (1. - fV(tkt, rain)) * taur) / rhow &
         + (fV(0., rain) - fV(tkt, rain)) * usrk**2)
    usrcs = sqrt((tau + (1. - fV(tks, rain)) * taur) / rhow &
         + (fV(0., rain) - fV(tks, rain)) * usrk**2)

    where (alq > 0.)
       ! Fairall 1996 982, equation (14):
       xlamx = 6. * (1. + (16. * grav * cpw * rhow * visw**3 * alq &
            / (tcw**2 * usrct**4 ))**0.75)**(- 1. / 3.)

       ! Fairall 1996 982, equation (12):
       tkt = xlamx * visw / usrct

       tks = xlamx * mu * (kappa / mu)**(2. / 3.) &
            * visw * cpw * rhow / ( tcw * usrcs)
       ! From Saunders 1967 (4)
    elsewhere
       xlamx = 6. ! prevent excessive warm skins
       tkt = min(.01, xlamx * visw / usrct) ! Limit tkt
       tks = min(.001, xlamx * mu * (kappa / mu)**(2. / 3.) * visw * cpw &
            * rhow / (tcw * usrcs))
    end where

    ! Fairall 1996 982, equation (13):
    dter = - (qcol + rf * (1 - fV(tkt, rain))) * tkt / tcw

    dser = s_subskin * (hlb / xlv - rain * (1 - fV(tks, rain))) * tks &
         / (rhow * mu) ! eq. fresh skin

  end subroutine Microlayer

end module Microlayer_m