source: LMDZ6/branches/Amaury_dev/libf/phylmd/Ocean_skin/near_surface_m.F90

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