module bulk_flux_m

  implicit none

contains

  subroutine bulk_flux(tkt, tks, taur, dter, dser, t_int, s_int, ds_ns, dt_ns, &
       u, t_ocean_1, s1, rain, hf, hlb, rnl, tau, rhoa, xlv, rf, dtime, rns)

    use config_ocean_skin_m, only: jwarm, jcool, rain_effect
    use Microlayer_m, only: Microlayer
    use mom_flux_rain_m, only: mom_flux_rain
    use Near_Surface_m, only: Near_Surface, depth
    use therm_expans_m, only: therm_expans

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

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

    real, intent(out):: taur(:) ! momentum flux due to rain, in Pa

    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(out):: t_int(:) ! interface temperature, in K
    real, intent(out):: s_int(:) ! interface salinity, in ppt

    real, intent(inout):: ds_ns(:)
    ! "delta salinity near surface". Salinity variation in the
    ! near-surface turbulent layer. That is subskin salinity minus
    ! foundation salinity. In ppt.

    real, intent(inout):: dt_ns(:)
    ! "delta temperature near surface". Temperature variation in the
    ! near-surface turbulent layer. That is subskin temperature minus
    ! foundation temperature. (Can be negative.) In K.

    real, intent(in):: u(:)
    ! Wind speed relative to the sea surface, i. e. taking current
    ! vector into account. In m s-1.

    real, intent(in):: t_ocean_1(:) ! input sea temperature, at depth_1, in K
    real, intent(in):: S1(:) ! salinity at depth_1, in ppt

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

    real, intent(in):: hf(:)
    !  turbulent part of sensible heat flux, positive upward, in W m-2

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

    real, intent(in):: rnl(:)
    ! net longwave radiation, positive upward, in W m-2

    real, intent(in):: tau(:)
    ! wind stress at the surface, turbulent part only, in Pa
    
    real, intent(in):: rhoa(:) ! density of moist air  (kg / m3)
    real, intent(in):: xlv(:) ! latent heat of evaporation (J / kg)

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

    real, intent(in):: dtime ! time step, in s
    real, intent(in):: rns(:) ! net downward shortwave radiation, in W m-2

    ! Local:
    
    real al(size(t_ocean_1)) ! water thermal expansion coefficient (in K-1)
    real dels(size(t_ocean_1)), null_array(size(t_ocean_1))
    integer iter
    real t_subskin(size(t_ocean_1)) ! subskin temperature, in K
    real s_subskin(size(t_ocean_1)) ! subskin salinity, in ppt

    real, parameter:: fxp = 1. - (0.28 * 0.014 &
         + 0.27 * 0.357 * (1. - exp(- depth / 0.357)) &
         + .45 * 12.82 * (1.- exp(- depth / 12.82))) / depth
    ! underflow ! fxp = 1. - (0.28 * 0.014 * (1. - exp(- depth / 0.014)) &
    ! Soloviev solar absorption profile
    ! H. Bellenger 2016

    real tau_with_min(size(t_ocean_1))
    ! modified wind stress, avoiding very low values
    
    real, parameter:: tau_0 = 1e-3 ! in N m-2

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

    if (rain_effect) then
       taur = mom_flux_rain(u, rain)
    else
       if (jwarm .or. jcool) null_array = 0.
       taur = 0.
    end if

    if (jwarm .or. jcool) tau_with_min = tau + tau_0 * (1. - exp(- tau_0 / tau))

    if (Jwarm) then
       if (rain_effect) then
          call Near_Surface(al, t_subskin, s_subskin, ds_ns, dt_ns, &
               tau_with_min, taur, hlb, rhoa, xlv, dtime, t_ocean_1, s1, rain, &
               q_pwp = fxp * rns - (hf + hlb + rnl + rf))
       else
          call Near_Surface(al, t_subskin, s_subskin, ds_ns, dt_ns, &
               tau_with_min, taur, hlb, rhoa, xlv, dtime, t_ocean_1, s1, &
               rain = null_array, q_pwp = fxp * rns - (hf + hlb + rnl))
       end if
    else
       if (Jcool) al = therm_expans(t_ocean_1)
       t_subskin = t_ocean_1
       s_subskin = s1
    end if

    if (Jcool) then
       ! First guess:
       tkt = 0.001
       tks = 5e-4

       do iter = 1, 3
          ! Cool skin
          dels = rns * (0.065 + 11. * tkt - 6.6e-5 / tkt &
               * (1. - exp(- tkt / 8e-4))) ! equation 16 Ohlmann
          if (rain_effect) then
             call Microlayer(dter, dser, tkt, tks, hlb, tau_with_min, &
                  s_subskin, al, xlv, taur, rf, rain, &
                  qcol = rnl + hf + hlb - dels)
          else
             call Microlayer(dter, dser, tkt, tks, hlb, tau_with_min, &
                  s_subskin, al, xlv, taur, rf = null_array, &
                  rain = null_array, qcol = rnl + hf + hlb - dels)
          end if
       end do
    else
       tkt = 0.
       tks = 0.
       dter = 0.
       dser = 0.
    end if

    t_int = t_subskin + dter
    s_int = s_subskin + dser

  end subroutine bulk_flux

end module bulk_flux_m