MODULE lmdz_spla_nightingale
CONTAINS
SUBROUTINE spla_nightingale(klon,klev,nbsrf,u, v, u_10m, v_10m, paprs, pplay, &
        cdragh, cdragm, t, q, ftsol, tsol, &
        pctsrf, lmt_dmsconc, lmt_dms)
  !
USE lmdz_spla_ini, ONLY: is_oce, RNAVO

IMPLICIT NONE
  !
INTEGER klon,klev,nbsrf
  !
  REAL, dimension(klon,klev), intent(in) :: u, v
  REAL, dimension(klon), intent(in):: u_10m, v_10m
  REAL, dimension(klon,nbsrf), intent(in):: ftsol
  REAL, dimension(klon), intent(in) :: tsol
  REAL, dimension(klon,klev+1), intent(in) :: paprs
  REAL, dimension(klon,klev), intent(in) :: pplay
  REAL, dimension(klon,klev), intent(in) :: t
  REAL, dimension(klon,klev), intent(in) :: q
  REAL, dimension(klon), intent(in) :: cdragh, cdragm
  REAL, dimension(klon,nbsrf), intent(in) :: pctsrf
  REAL, dimension(klon), intent(out) :: lmt_dmsconc  ! concentration oceanique DMS
  REAL, dimension(klon), intent(out) :: lmt_dms      ! flux de DMS
  !
  REAL, dimension(klon) :: ustar, obklen
  REAL, dimension(klon) :: u10, u10n
  REAL :: tvelocity, schmidt_corr
  REAL :: t1, t2, t3, t4, viscosity_kin, diffusivity, schmidt
  INTEGER :: i
  !
  CALL bl_for_dms(u, v, paprs, pplay, cdragh, cdragm, &
        t, q, tsol, ustar, obklen)
  !
  DO i=1,klon
    u10(i)=SQRT(u_10m(i)**2+v_10m(i)**2)
  ENDDO
  !
  CALL neutral(u10, ustar, obklen, u10n)
  !
  DO i=1,klon
  !
  ! tvelocity - transfer velocity, also known as kw (cm/s)
  ! schmidt_corr - Schmidt number correction factor (dimensionless)
  ! Reference:  Nightingale, P.D., G. Malin, C. S. Law, J. J. Watson, P.S. Liss
  !  M. I. Liddicoat, J. Boutin, R.C. Upstill-Goddard. 'In situ evaluation
  !  of air-sea gas exchange parameterizations using conservative and
  !  volatile tracers.'  Glob. Biogeochem. Cycles, 14:373-387, 2000.
  ! compute transfer velocity using u10neutral
  !
  tvelocity = 0.222*u10n(i)*u10n(i) + 0.333*u10n(i)
  !
  ! above expression gives tvelocity in cm/hr. convert to cm/s. 1hr =3600 sec

  tvelocity = tvelocity / 3600.

  ! compute the correction factor, which for Nightingale parameterization is
  ! based on how different the schmidt number is from 600.
  ! correction factor based on temperature in Kelvin. good
  ! only for t<=30 deg C.  for temperatures above that, set correction factor
  ! equal to value at 30 deg C.

  IF (ftsol(i,is_oce) .LE. 303.15) THEN
     t1 = ftsol(i,is_oce)
  ELSE
     t1 = 303.15
  ENDIF

  t2 = t1 * t1
  t3 = t2 * t1
  t4 = t3 * t1
  viscosity_kin = 3.0363e-9*t4 - 3.655198e-6*t3 + 1.65333e-3*t2 &
        - 3.332083e-1*t1 + 25.26819
  diffusivity = 0.01922 * exp(-2177.1/t1)
  schmidt = viscosity_kin / diffusivity
  schmidt_corr = (schmidt/600.)**(-.5)
  !
  lmt_dms(i) = tvelocity  *  pctsrf(i,is_oce) &
        * lmt_dmsconc(i)/1.0e12 * schmidt_corr * RNAVO
  !
  IF (lmt_dmsconc(i).LE.1.e-20) lmt_dms(i)=0.0
  !
  ENDDO
  !
END SUBROUTINE spla_nightingale
END MODULE lmdz_spla_nightingale