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

    ! 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.

    ! Wind speed at 15 m relative to the sea surface, i. e. taking

    ! current vector into account. In m s-1.

    ! exponent in the function giving T(z) and S(z), equation (11) in

    ! Bellenger et al. 2017 JGR

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

    ! Temperature and salinity profiles change with wind:

       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

       ! (neglecting microlayer thickness compared to depth_1 and depth)

    else

    end if

    ! Bellenger 2017 k0976, equation (13):

    buoyf = Al * grav / (rhow * cpw) * q_pwp &

    ! Case of stable stratification and negative flux, Bellenger 2017

    ! k0976, equation (15):

       buoyf = sqrt(- eta * grav / (5. * depth * rhow) * drho) * usrc**2

    elsewhere (buoyf == 0.)

       buoyf = tiny(0.)

    end where

    ! Equation (14) for temperature. Implicit scheme for time integration:

    ! \Delta T_{i + 1} - \Delta T_i = \delta t (Bt + At \Delta T_{i + 1})

    ! Lens horizontal spreading:

    where (drho < 0. .and. ds_ns < 0.) At = At &

    ! Equation (14) for salinity:

    ! \frac{\partial \Delta S}{\partial t}

    ! = (\Delta S + S_\mathrm{fnd}) B_S + A_S \Delta S

    ! Lens horizontal spreading:

    where (drho < 0. .and. ds_ns < 0.) As = As &

    ! Implicit scheme for time integration:

end module Near_Surface_m