MODULE lmdz_spla_bl_for_dms

CONTAINS

SUBROUTINE spla_bl_for_dms(klon,klev,u,v,paprs,pplay,cdragh,cdragm &
        ,t,q,tsol,ustar,obklen)

  USE lmdz_spla_ini, ONLY : RG,RD,RCPD,RLVTT,RLSTT,RETV,RTT, &
          R2ES,R3IES,R3LES,R4IES,R4LES,R5IES,R5LES,RVTMP2

 !USE yoethf_mod_h
  !
  !===================================================================
  ! Auteur : E. Cosme
  ! Calcul de la vitesse de friction (ustar) et de la longueur de
  ! Monin-Obukhov (obklen), necessaires pour calculer les flux de DMS
  ! par la methode de Nightingale.
  ! Cette subroutine est plus que fortement inspiree de la subroutine
  ! 'nonlocal' dans clmain.F .
  ! reference :  Holtslag, A.A.M., and B.A. Boville, 1993:
  ! Local versus nonlocal boundary-layer diffusion in a global climate
  ! model. J. of Climate, vol. 6, 1825-1842. (a confirmer)
  ! 31 08 01
  !===================================================================
  !
  IMPLICIT NONE

  INCLUDE "FCTTRE.h"

  !
  INTEGER, intent(in) :: klon,klev
  !
  ! Arguments :
  REAL, dimension(klon,klev), intent(in) :: u          ! vent zonal
  REAL, dimension(klon,klev), intent(in) :: v          ! vent meridien
  REAL, dimension(klon,klev+1), intent(in):: paprs    ! niveaux de pression aux intercouches (Pa)
  REAL, dimension(klon,klev), intent(in) :: pplay     ! niveaux de pression aux milieux... (Pa)
  REAL, dimension(klon,klev), intent(in) :: t          ! temperature
  REAL, dimension(klon,klev), intent(in) :: q          ! humidite kg/kg
  REAL, dimension(klon), intent(in) :: cdragh          ! coefficient de trainee pour la chaleur
  REAL, dimension(klon), intent(in) :: cdragm          ! coefficient de trainee pour le vent
  REAL, dimension(klon), intent(in) :: tsol            ! temperature du sol
  REAL, dimension(klon), intent(out) :: ustar           ! vitesse de friction
  REAL, dimension(klon), intent(out) :: obklen          ! longueur de Monin-Obukhov
  !
  ! Locales :
  REAL :: vk
  PARAMETER (vk=0.35)
  REAL :: beta  ! coefficient d'evaporation reelle (/evapotranspiration)
             ! ! entre 0 et 1, mais 1 au-dessus de la mer
  PARAMETER (beta=1.)
  INTEGER :: i,k
  REAL :: zxt, zxu, zxv, zxq, zxqs, zxmod, taux, tauy
  REAL :: zcor, zdelta, zcvm5
  REAL, dimension(klon,klev) :: z
  REAL :: zx_alf1, zx_alf2 ! parametres pour extrapolation
  REAL, dimension(klon) :: khfs       ! surface kinematic heat flux [mK/s]
  REAL, dimension(klon) :: kqfs       ! sfc kinematic constituent flux [m/s]
  REAL, dimension(klon) :: heatv      ! surface virtual heat flux


  !
  !======================================================================
  !
  ! Calculer les hauteurs de chaque couche
  !
  ! JE20150707      r2es=611.14 *18.0153/28.9644
  DO i = 1, klon
     z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
           * (paprs(i,1)-pplay(i,1)) / RG
  ENDDO
  DO k = 2, klev
  DO i = 1, klon
     z(i,k) = z(i,k-1) &
           + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) &
           * (pplay(i,k-1)-pplay(i,k)) / RG
  ENDDO
  ENDDO

  DO i = 1, klon
  !
    zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
    zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
    zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
    zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
    zxqs=MIN(0.5,zxqs)
    zcor=1./(1.-retv*zxqs)
    zxqs=zxqs*zcor
  !
    zx_alf1 = 1.0
    zx_alf2 = 1.0 - zx_alf1
    zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1))) &
          *(1.+RETV*q(i,1))*zx_alf1 &
          + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2))) &
          *(1.+RETV*q(i,2))*zx_alf2
    zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
    zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
    zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
    zxmod = 1.0+SQRT(zxu**2+zxv**2)
    khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cdragh(i)
    kqfs(i) = (zxqs-zxq) *zxmod*cdragh(i) * beta
    heatv(i) = khfs(i) + 0.61*zxt*kqfs(i)
    taux = zxu *zxmod*cdragm(i)
    tauy = zxv *zxmod*cdragm(i)
    ustar(i) = SQRT(taux**2+tauy**2)
    ustar(i) = MAX(SQRT(ustar(i)),0.01)
  !
  ENDDO
  !
  DO i = 1, klon
     obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
  ENDDO
  !
END SUBROUTINE spla_bl_for_dms

END MODULE lmdz_spla_bl_for_dms