module Near_Surface_m

  Implicit none

  real, parameter:: depth = 3.
  ! diurnal warm layer and fresh water lens depth, in m (Zeng and Beljaars 2005)

contains

  subroutine near_surface(al, t_subskin, s_subskin, ds_ns, dt_ns, tau, taur, &
       hlb, rhoa, xlv, dtime, t_ocean_1, s1, rain, q_pwp)

    ! Hugo Bellenger, 2016

    use config_ocean_skin_m, only: depth_1
    use const, only: beta, cpw, grav, rhow, von
    use Phiw_m, only: Phiw
    use therm_expans_m, only: therm_expans

    real, intent(out):: al(:) ! water thermal expansion coefficient (in K-1)
    real, intent(out):: t_subskin(:) ! subskin temperature, in K
    real, intent(out):: s_subskin(:) ! subskin 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):: tau(:)
    ! wind stress at the surface, turbulent part only, in Pa

    real, intent(in):: taur(:) ! momentum flux due to rainfall, in Pa
    real, intent(in):: hlb(:) ! latent heat flux, turbulent part only, in W / m2
    real, intent(in):: rhoa(:) ! density of moist air  (kg / m3)
    real, intent(in):: xlv(:) ! latent heat of evaporation (J/kg)
    real, intent(in):: dtime ! time step (s)
    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, in kg m-2 s-1

    real, intent(in):: q_pwp(:)
    ! net flux absorbed by the warm layer (part of the solar flux
    ! absorbed at "depth"), minus surface fluxes, in W m-2

    ! Local:

    real, parameter:: khor = 1. / 1.5e4
    ! Parameter for the lens spread, in m-1. Inverse of the size of
    ! the lens.

    real, parameter:: umax = 15.
    real, parameter:: fact = 1.
    real buoyf(size(t_ocean_1)) ! buoyancy flux
    real usrc(size(t_ocean_1))
    real drho(size(t_ocean_1)) ! rho(- delta) - rho(- d)
    real Lmo(size(t_ocean_1)) ! Monin-Obukhov length

    real u(size(t_ocean_1))
    ! Wind speed at 15 m relative to the sea surface, i. e. taking
    ! current vector into account. In m s-1.

    real, dimension(size(t_ocean_1)):: At, Bt, As, Bs, correction

    real eta(size(t_ocean_1))
    ! exponent in the function giving T(z) and S(z), equation (11) in
    ! Bellenger et al. 2017 JGR

    real t_fnd(size(t_ocean_1)) ! foundation temperature, in K
    real s_fnd(size(t_ocean_1)) ! foundation salinity, in ppt

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

    ! Temperature and salinity profiles change with wind:

    u = 28. * sqrt(tau / rhoa)

    where (dt_ns < 0.)
       where (u >= umax)
          eta = 1. / fact
       elsewhere (u <= 2.)
          eta = 2. / (fact * umax)
       elsewhere
          ! {u > 2. .and. u < umax}
          eta = u / (fact * umax)
       end where
    elsewhere
       eta = 0.3
    end where

    if (depth_1 < depth) then
       correction = 1. - (depth_1 / depth)**eta
       ! (neglecting microlayer thickness compared to depth_1 and depth)

       t_fnd = t_ocean_1 - dt_ns * correction
       s_fnd = s1 - ds_ns * correction
    else
       t_fnd = t_ocean_1
       s_fnd = s1
    end if

    al = therm_expans(t_fnd)

    ! Bellenger 2017 k0976, equation (13):
    buoyf = Al * grav / (rhow * cpw) * q_pwp &
         - beta * S_FND * grav * (hlb / xlv - rain) / rhow

    usrc = sqrt((tau + taur) / rhow)
    drho = rhow * (- al * dt_ns + beta * ds_ns)

    ! Case of stable stratification and negative flux, Bellenger 2017
    ! k0976, equation (15):
    where (buoyf < 0. .and. drho < 0.)
       buoyf = sqrt(- eta * grav / (5. * depth * rhow) * drho) * usrc**2
    elsewhere (buoyf == 0.)
       buoyf = tiny(0.)
    end where
       

    Lmo = usrc**3 / (von * buoyf)

    ! Equation (14) for temperature. Implicit scheme for time integration:
    ! \Delta T_{i + 1} - \Delta T_i = \delta t (Bt + At \Delta T_{i + 1})
    At = - (eta + 1.) * von * usrc / (depth * Phiw(depth / Lmo))

    ! Lens horizontal spreading:
    where (drho < 0. .and. ds_ns < 0.) At = At &
         - (eta + 1.) * khor * sqrt(depth * grav * abs(drho) / rhow)

    Bt = q_pwp / (depth * rhow * cpw * eta / (eta + 1.))
    dt_ns = (dtime * Bt + dt_ns) / (1 - dtime * At)

    ! Equation (14) for salinity:
    ! \frac{\partial \Delta S}{\partial t}
    ! = (\Delta S + S_\mathrm{fnd}) B_S + A_S \Delta S
    As = - (eta + 1.) * von * usrc / (depth * Phiw(depth / Lmo))

    ! Lens horizontal spreading:
    where (drho < 0. .and. ds_ns < 0.) As = As &
         - (eta + 1.) * khor * sqrt(depth * grav * abs(drho) / rhow)

    Bs = (hlb / xlv - rain) * (eta + 1.) / (depth * rhow * eta)

    ! Implicit scheme for time integration:
    ds_ns = (dtime * Bs * S_fnd + ds_ns) / (1 - dtime * (As + bs))

    t_subskin = t_fnd + dt_ns
    s_subskin = s_fnd + ds_ns

  end subroutine Near_Surface

end module Near_Surface_m