Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/bulk_flux.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/bulk_flux.F90	(revision 3833)
+++ 	(revision )
@@ -1,155 +1,0 @@
-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
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/bulk_flux_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/bulk_flux_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/bulk_flux_m.F90	(revision 3834)
@@ -0,0 +1,155 @@
+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
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/config_ocean_skin.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/config_ocean_skin.F90	(revision 3833)
+++ 	(revision )
@@ -1,81 +1,0 @@
-module config_ocean_skin_m
-
-  implicit none
-
-  logical, protected:: jcool ! cool skin calculation
-  logical, protected:: jwarm ! warm layer calculation
-  logical, protected:: rain_effect
-
-  real, protected:: depth_1 = 20.
-  ! Depth at which the temperature and salinity are input. Could be
-  ! the depth of a sensor, or the depth at the middle of the first
-  ! layer of an ocean model (half the depth of the first layer). In
-  ! m. Setting depth_1 to any value >= depth has the same effect as
-  ! setting depth_1 to depth.
-
-  !$omp threadprivate(jcool, jwarm, rain_effect, depth_1)
-
-#ifdef IN_LMDZ
-  integer, protected:: activate_ocean_skin = 0
-  ! Allowed values: 0 do not activate; 1 activate without retroaction
-  ! on LMDZ; 2 activate with retroaction on LMDZ
-  
-  !$omp threadprivate(activate_ocean_skin)
-#endif
-
-contains
-
-  subroutine config_ocean_skin
-
-#ifdef IN_LMDZ
-    use ioipsl_getin_p_mod, only: getin_p
-    use assert_m, only: assert
-#endif
-
-    integer:: flag_ocean_skin = 3
-    !$omp threadprivate(flag_ocean_skin)
-    
-    namelist /config_ocean_skin_nml/ flag_ocean_skin, depth_1
-
-    !-----------------------------------------------------------------------
-
-#ifdef IN_LMDZ
-    call getin_p("activate_ocean_skin", activate_ocean_skin)
-    call assert(activate_ocean_skin >= 0 .and. activate_ocean_skin <= 2, &
-         "config_ocean_skin bad value of activate_ocean_skin")
-    if (activate_ocean_skin >= 1) then
-       call getin_p("flag_ocean_skin", flag_ocean_skin)
-       call getin_p("depth_1", depth_1)
-    end if
-#else
-    write(unit = *, nml = config_ocean_skin_nml)
-    print *, "Enter namelist config_ocean_skin_nml."
-    read(unit = *, nml = config_ocean_skin_nml)
-    write(unit = *, nml = config_ocean_skin_nml)
-#endif
-
-    select case (flag_ocean_skin)
-    case (0)
-       jwarm = .false.
-       jcool = .false.
-       rain_effect = .false.
-    case (1)
-       jwarm = .false.
-       jcool = .true.
-       rain_effect = .false.
-    case (2)
-       jwarm = .true.
-       jcool = .true.
-       rain_effect = .false.
-    case (3)
-       jwarm = .true.
-       jcool = .true.
-       rain_effect = .true.
-    case default
-       print *, "config_ocean_skin: bad value for flag_ocean_skin."
-       stop 1
-    end select
-
-  end subroutine config_ocean_skin
-
-end module config_ocean_skin_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/config_ocean_skin_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/config_ocean_skin_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/config_ocean_skin_m.F90	(revision 3834)
@@ -0,0 +1,81 @@
+module config_ocean_skin_m
+
+  implicit none
+
+  logical, protected:: jcool ! cool skin calculation
+  logical, protected:: jwarm ! warm layer calculation
+  logical, protected:: rain_effect
+
+  real, protected:: depth_1 = 20.
+  ! Depth at which the temperature and salinity are input. Could be
+  ! the depth of a sensor, or the depth at the middle of the first
+  ! layer of an ocean model (half the depth of the first layer). In
+  ! m. Setting depth_1 to any value >= depth has the same effect as
+  ! setting depth_1 to depth.
+
+  !$omp threadprivate(jcool, jwarm, rain_effect, depth_1)
+
+#ifdef IN_LMDZ
+  integer, protected:: activate_ocean_skin = 0
+  ! Allowed values: 0 do not activate; 1 activate without retroaction
+  ! on LMDZ; 2 activate with retroaction on LMDZ
+  
+  !$omp threadprivate(activate_ocean_skin)
+#endif
+
+contains
+
+  subroutine config_ocean_skin
+
+#ifdef IN_LMDZ
+    use ioipsl_getin_p_mod, only: getin_p
+    use assert_m, only: assert
+#endif
+
+    integer:: flag_ocean_skin = 3
+    !$omp threadprivate(flag_ocean_skin)
+    
+    namelist /config_ocean_skin_nml/ flag_ocean_skin, depth_1
+
+    !-----------------------------------------------------------------------
+
+#ifdef IN_LMDZ
+    call getin_p("activate_ocean_skin", activate_ocean_skin)
+    call assert(activate_ocean_skin >= 0 .and. activate_ocean_skin <= 2, &
+         "config_ocean_skin bad value of activate_ocean_skin")
+    if (activate_ocean_skin >= 1) then
+       call getin_p("flag_ocean_skin", flag_ocean_skin)
+       call getin_p("depth_1", depth_1)
+    end if
+#else
+    write(unit = *, nml = config_ocean_skin_nml)
+    print *, "Enter namelist config_ocean_skin_nml."
+    read(unit = *, nml = config_ocean_skin_nml)
+    write(unit = *, nml = config_ocean_skin_nml)
+#endif
+
+    select case (flag_ocean_skin)
+    case (0)
+       jwarm = .false.
+       jcool = .false.
+       rain_effect = .false.
+    case (1)
+       jwarm = .false.
+       jcool = .true.
+       rain_effect = .false.
+    case (2)
+       jwarm = .true.
+       jcool = .true.
+       rain_effect = .false.
+    case (3)
+       jwarm = .true.
+       jcool = .true.
+       rain_effect = .true.
+    case default
+       print *, "config_ocean_skin: bad value for flag_ocean_skin."
+       stop 1
+    end select
+
+  end subroutine config_ocean_skin
+
+end module config_ocean_skin_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/esat.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/esat.F90	(revision 3833)
+++ 	(revision )
@@ -1,22 +1,0 @@
-module esat_m
-
-  implicit none
-
-contains
-
-  elemental real function esat(T, P)
-
-    ! Saturation vapor pressure of water in Pa. Buck, 1981,
-    ! J. Appl. Meteor. 20, 1527-1532, equation (8).
-
-    real, intent(in):: T ! temperature, in K
-    real, intent(in):: P ! air pressure, in Pa
-
-    !--------------------------------------------------------------------
-
-    esat = (1.0007 + 3.46e-8 * P) * 6.1121e2 &
-         * exp(17.502 * (T - 273.15) / (T - 32.18))
-
-  end function esat
-
-end module esat_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/esat_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/esat_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/esat_m.F90	(revision 3834)
@@ -0,0 +1,22 @@
+module esat_m
+
+  implicit none
+
+contains
+
+  elemental real function esat(T, P)
+
+    ! Saturation vapor pressure of water in Pa. Buck, 1981,
+    ! J. Appl. Meteor. 20, 1527-1532, equation (8).
+
+    real, intent(in):: T ! temperature, in K
+    real, intent(in):: P ! air pressure, in Pa
+
+    !--------------------------------------------------------------------
+
+    esat = (1.0007 + 3.46e-8 * P) * 6.1121e2 &
+         * exp(17.502 * (T - 273.15) / (T - 32.18))
+
+  end function esat
+
+end module esat_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/fv.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/fv.F90	(revision 3833)
+++ 	(revision )
@@ -1,46 +1,0 @@
-module fv_m
-
-  implicit none
-
-contains
-
-  elemental real function fV(z, rain)
-
-    ! H. Bellenger 2016
-
-    real, intent(in):: z ! profondeur en m (z < 0)
-    real, intent(in):: rain ! rain mass flux, in kg m-2 s-1
-
-    ! Local:
-
-    real dfV
-    ! fraction of rain volume entering the ocean and deposited within depth "z"
-
-    real rc, lbd, f0
-    real z_mm
-
-    ! Schlussel et al. 1997, Table 1:
-    real, parameter:: zeta(14) = [0.0026091, 0.0022743, 0.0015406, 0.0012281, &
-         0.0008795, 0.00077123, 0.00057451, 0.000438, 6.7228e-5, 6.4955e-5, &
-         4.4234e-5, 3.3906e-5, 2.7433e-6, 4.0283e-7]
-    real, parameter:: psi(14) = [0.60107, 0.29968, 0.5563, 1.80858, 0.2175, &
-         0.33961, 0.96368, 0.65081, 0.5967, 2.7661, 2.2812, 2.7674, 2.6095, &
-         6.5308]
-
-    !---------------------------------------------------------------------
-
-    if (rain > 0.) then
-       rc = 0.4 ! mm
-       z_mm = z * 1000. ! mm and <0
-       lbd = 4.1 * (rain * 3600.)**(- 0.21) ! mm-1
-       f0 = (1. + 2. * lbd * rc + 0.5 * (2. * lbd * rc)**2 &
-            + 1. / 6. * (2. * lbd * rc)**3) * exp(- 1. * 2. * lbd * rc)
-       dfv = sum(zeta * exp(- psi * lbd * abs(z_mm) / 100.))
-       fV = - dfV * abs(z_mm) + f0
-    else
-       fV = 0.
-    endif
-
-  end function fV
-
-end module fv_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/fv_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/fv_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/fv_m.F90	(revision 3834)
@@ -0,0 +1,46 @@
+module fv_m
+
+  implicit none
+
+contains
+
+  elemental real function fV(z, rain)
+
+    ! H. Bellenger 2016
+
+    real, intent(in):: z ! profondeur en m (z < 0)
+    real, intent(in):: rain ! rain mass flux, in kg m-2 s-1
+
+    ! Local:
+
+    real dfV
+    ! fraction of rain volume entering the ocean and deposited within depth "z"
+
+    real rc, lbd, f0
+    real z_mm
+
+    ! Schlussel et al. 1997, Table 1:
+    real, parameter:: zeta(14) = [0.0026091, 0.0022743, 0.0015406, 0.0012281, &
+         0.0008795, 0.00077123, 0.00057451, 0.000438, 6.7228e-5, 6.4955e-5, &
+         4.4234e-5, 3.3906e-5, 2.7433e-6, 4.0283e-7]
+    real, parameter:: psi(14) = [0.60107, 0.29968, 0.5563, 1.80858, 0.2175, &
+         0.33961, 0.96368, 0.65081, 0.5967, 2.7661, 2.2812, 2.7674, 2.6095, &
+         6.5308]
+
+    !---------------------------------------------------------------------
+
+    if (rain > 0.) then
+       rc = 0.4 ! mm
+       z_mm = z * 1000. ! mm and <0
+       lbd = 4.1 * (rain * 3600.)**(- 0.21) ! mm-1
+       f0 = (1. + 2. * lbd * rc + 0.5 * (2. * lbd * rc)**2 &
+            + 1. / 6. * (2. * lbd * rc)**3) * exp(- 1. * 2. * lbd * rc)
+       dfv = sum(zeta * exp(- psi * lbd * abs(z_mm) / 100.))
+       fV = - dfV * abs(z_mm) + f0
+    else
+       fV = 0.
+    endif
+
+  end function fV
+
+end module fv_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/microlayer.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/microlayer.F90	(revision 3833)
+++ 	(revision )
@@ -1,106 +1,0 @@
-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
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/microlayer_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/microlayer_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/microlayer_m.F90	(revision 3834)
@@ -0,0 +1,106 @@
+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
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/mom_flux_rain.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/mom_flux_rain.F90	(revision 3833)
+++ 	(revision )
@@ -1,24 +1,0 @@
-module mom_flux_rain_m
-
-  implicit none
-
-contains
-
-  elemental real function mom_flux_rain(u, rain)
-
-    ! Computes momentum flux due to rainfall, in Pa.
-
-    real, intent(in):: U
-    ! difference of velocity between air and sea, including gustiness, in m / s
-
-    real, intent(in):: rain ! rain mass flux, in kg m-2 s-1
-
-    !------------------------------------------------------------------------
-
-    mom_flux_rain = 0.85 * U * rain
-    ! (Caldwell 1971 k1001, equation (1) and 15 % reduction in speed
-    ! as in paragraph 3.a, maybe adequate if u is the wind at 10 m)
-
-  end function mom_flux_rain
-  
-end module mom_flux_rain_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/mom_flux_rain_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/mom_flux_rain_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/mom_flux_rain_m.F90	(revision 3834)
@@ -0,0 +1,24 @@
+module mom_flux_rain_m
+
+  implicit none
+
+contains
+
+  elemental real function mom_flux_rain(u, rain)
+
+    ! Computes momentum flux due to rainfall, in Pa.
+
+    real, intent(in):: U
+    ! difference of velocity between air and sea, including gustiness, in m / s
+
+    real, intent(in):: rain ! rain mass flux, in kg m-2 s-1
+
+    !------------------------------------------------------------------------
+
+    mom_flux_rain = 0.85 * U * rain
+    ! (Caldwell 1971 k1001, equation (1) and 15 % reduction in speed
+    ! as in paragraph 3.a, maybe adequate if u is the wind at 10 m)
+
+  end function mom_flux_rain
+  
+end module mom_flux_rain_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/near_surface.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/near_surface.F90	(revision 3833)
+++ 	(revision )
@@ -1,156 +1,0 @@
-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
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/near_surface_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/near_surface_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/near_surface_m.F90	(revision 3834)
@@ -0,0 +1,156 @@
+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
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/phiw.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/phiw.F90	(revision 3833)
+++ 	(revision )
@@ -1,24 +1,0 @@
-module Phiw_m
-
-  implicit none
-
-contains
-
-  elemental real function Phiw(zL)
-
-    ! Fonction de stabilit\'e (Monin Obukhov)
-
-    real, intent(in):: zL ! profondeur  / longueur de Monin-Obukhov
-
-    !----------------------------------------------
-
-    if (zL < 0.) then
-       Phiw = (1. - 16. * zL)**(- 0.5)
-    else
-       Phiw = 1. + (5. * zL + 4. * zL**2)/(1. + 3. * zL + 0.25 * zL**2)
-       ! (Takaya et al. 2010)
-    end if
-
-  end function Phiw
-
-end module Phiw_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/phiw_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/phiw_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/phiw_m.F90	(revision 3834)
@@ -0,0 +1,24 @@
+module Phiw_m
+
+  implicit none
+
+contains
+
+  elemental real function Phiw(zL)
+
+    ! Fonction de stabilit\'e (Monin Obukhov)
+
+    real, intent(in):: zL ! profondeur  / longueur de Monin-Obukhov
+
+    !----------------------------------------------
+
+    if (zL < 0.) then
+       Phiw = (1. - 16. * zL)**(- 0.5)
+    else
+       Phiw = 1. + (5. * zL + 4. * zL**2)/(1. + 3. * zL + 0.25 * zL**2)
+       ! (Takaya et al. 2010)
+    end if
+
+  end function Phiw
+
+end module Phiw_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/sens_heat_rain.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/sens_heat_rain.F90	(revision 3833)
+++ 	(revision )
@@ -1,66 +1,0 @@
-module sens_heat_rain_m
-
-  implicit none
-
-contains
-
-  real function sens_heat_rain(rain, t, q, rhoa, xlv, t_int, p)
-
-    ! Computes heat flux due to rainfall, in W m-2, positive
-    ! upward.
-
-    ! Do not declare this function as elemental because it may include
-    ! YOMCST.h, which contains an OpenMP directive, and that only
-    ! works in OpenMP 5.0.
-
-    use const, only: cpa, cpw, rgas
-#ifndef IN_LMDZ
-    use const, only: eps_w
-#endif
-    use esat_m, only: esat
-
-#ifdef IN_LMDZ
-    include "YOMCST.h"
-    ! for eps_w
-#endif
-
-    real, intent(in):: rain ! rain mass flux, in kg m-2 s-1
-    real, intent(in):: t ! air temperature at 10 m, in K
-    real, intent(in):: q ! specific humidity at 10 m
-    real, intent(in):: rhoa ! density of moist air  (kg / m3)
-    real, intent(in):: xlv ! latent heat of evaporation (J / kg)
-    real, intent(in):: t_int ! interface temperature, in K
-    real, intent(in):: p ! surface pressure, in Pa
-
-    ! Local:
-    
-    real es ! saturation pressure of wator vapor, in Pa
-    real alfac ! wet bulb factor
-    real dwat ! water vapour diffusivity
-    real dtmp ! heat diffusivity
-    real q_int ! specific (saturation) humidity at ocean interface
-    real t_celsius ! air temperature at 10 m, in Celsius degrees
-
-    real wetc
-    ! derivative of saturated mass fraction of water vapor with
-    ! respect to temperature, at constant total pressure
-
-    !---------------------------------------------------------------------
-
-    es = esat(t_int, p) * 0.98 ! reduced for salinity, Kraus 1972 page 46
-    q_int = eps_w * (es / (p - (1. - eps_w) * es))
-    wetc = eps_w * xlv * q_int / (rgas * t_int**2)
-    dwat = 2.11e-5 * (t / 273.15)**1.94
-    t_celsius = t - 273.15
-    dtmp = (1. + 3.309e-3 * t_celsius - 1.44e-6 * t_celsius**2) * 0.02411 &
-         / (rhoa * cpa)
-
-    ! Gosnell 1995 k0991, equation (11):
-    alfac =  1. / (1. + (wetc * xlv * dwat) / (cpa * dtmp))
-
-    ! Gosnell 1995 k0991, equation (12):
-    sens_heat_rain =  rain * alfac * cpw * (t_int - t + (q_int - q) * xlv / cpa)
-
-  end function sens_heat_rain
-  
-end module sens_heat_rain_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/sens_heat_rain_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/sens_heat_rain_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/sens_heat_rain_m.F90	(revision 3834)
@@ -0,0 +1,66 @@
+module sens_heat_rain_m
+
+  implicit none
+
+contains
+
+  real function sens_heat_rain(rain, t, q, rhoa, xlv, t_int, p)
+
+    ! Computes heat flux due to rainfall, in W m-2, positive
+    ! upward.
+
+    ! Do not declare this function as elemental because it may include
+    ! YOMCST.h, which contains an OpenMP directive, and that only
+    ! works in OpenMP 5.0.
+
+    use const, only: cpa, cpw, rgas
+#ifndef IN_LMDZ
+    use const, only: eps_w
+#endif
+    use esat_m, only: esat
+
+#ifdef IN_LMDZ
+    include "YOMCST.h"
+    ! for eps_w
+#endif
+
+    real, intent(in):: rain ! rain mass flux, in kg m-2 s-1
+    real, intent(in):: t ! air temperature at 10 m, in K
+    real, intent(in):: q ! specific humidity at 10 m
+    real, intent(in):: rhoa ! density of moist air  (kg / m3)
+    real, intent(in):: xlv ! latent heat of evaporation (J / kg)
+    real, intent(in):: t_int ! interface temperature, in K
+    real, intent(in):: p ! surface pressure, in Pa
+
+    ! Local:
+    
+    real es ! saturation pressure of wator vapor, in Pa
+    real alfac ! wet bulb factor
+    real dwat ! water vapour diffusivity
+    real dtmp ! heat diffusivity
+    real q_int ! specific (saturation) humidity at ocean interface
+    real t_celsius ! air temperature at 10 m, in Celsius degrees
+
+    real wetc
+    ! derivative of saturated mass fraction of water vapor with
+    ! respect to temperature, at constant total pressure
+
+    !---------------------------------------------------------------------
+
+    es = esat(t_int, p) * 0.98 ! reduced for salinity, Kraus 1972 page 46
+    q_int = eps_w * (es / (p - (1. - eps_w) * es))
+    wetc = eps_w * xlv * q_int / (rgas * t_int**2)
+    dwat = 2.11e-5 * (t / 273.15)**1.94
+    t_celsius = t - 273.15
+    dtmp = (1. + 3.309e-3 * t_celsius - 1.44e-6 * t_celsius**2) * 0.02411 &
+         / (rhoa * cpa)
+
+    ! Gosnell 1995 k0991, equation (11):
+    alfac =  1. / (1. + (wetc * xlv * dwat) / (cpa * dtmp))
+
+    ! Gosnell 1995 k0991, equation (12):
+    sens_heat_rain =  rain * alfac * cpw * (t_int - t + (q_int - q) * xlv / cpa)
+
+  end function sens_heat_rain
+  
+end module sens_heat_rain_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/therm_expans.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/therm_expans.F90	(revision 3833)
+++ 	(revision )
@@ -1,19 +1,0 @@
-module therm_expans_m
-
-  implicit none
-
-contains
-
-  elemental real function therm_expans(t)
-
-    ! Compute the thermal expansion coefficient of sea water.
-
-    real, intent(in):: t ! temperature, in K
-
-    !--------------------------------------------------------------
-
-    therm_expans = 2.1e-5 * (t - 269.95)**0.79
-
-  end function therm_expans
-
-end module therm_expans_m
Index: LMDZ6/trunk/libf/phylmd/Ocean_skin/therm_expans_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/Ocean_skin/therm_expans_m.F90	(revision 3834)
+++ LMDZ6/trunk/libf/phylmd/Ocean_skin/therm_expans_m.F90	(revision 3834)
@@ -0,0 +1,19 @@
+module therm_expans_m
+
+  implicit none
+
+contains
+
+  elemental real function therm_expans(t)
+
+    ! Compute the thermal expansion coefficient of sea water.
+
+    real, intent(in):: t ! temperature, in K
+
+    !--------------------------------------------------------------
+
+    therm_expans = 2.1e-5 * (t - 269.95)**0.79
+
+  end function therm_expans
+
+end module therm_expans_m
