MODULE lmdz_surf_wind
        CONTAINS 

SUBROUTINE surf_wind(klon,nsurfwind,zu10m,zv10m,sigmaw,cstar,ustar,ale_bl,wind10ms,probu)

USE lmdz_surf_wind_ini, ONLY : iflag_surf_wind
USE lmdz_surf_wind_ini, ONLY : surf_wind_ktwake
USE lmdz_surf_wind_ini, ONLY : surf_wind_kttherm
USE lmdz_surf_wind_ini, ONLY : surf_wind_kztherm

IMPLICIT NONE
INTEGER, INTENT(IN)                :: nsurfwind, klon
REAL, DIMENSION(klon), INTENT(IN)  :: zu10m, zv10m
REAL, DIMENSION(klon), INTENT(IN)  :: cstar
REAL, DIMENSION(klon), INTENT(IN)  :: sigmaw
REAL, DIMENSION(klon), INTENT(IN)  :: ustar, ale_bl
REAL, DIMENSION(klon,nsurfwind), INTENT(OUT)         :: wind10ms, probu
REAL, PARAMETER                    :: woff=0.5  ! min value of 10m wind speed accepted for emissions

REAL, DIMENSION(klon,nsurfwind)         :: sigma_th, sigma_wk
REAL, DIMENSION(klon,nsurfwind)         :: xp, yp
REAL, DIMENSION(klon,nsurfwind)         :: vwx, vwy, vw
REAL, DIMENSION(klon,nsurfwind)         :: vtx, vty
REAL, DIMENSION(klon,nsurfwind)         :: windx, windy, wind
REAL, DIMENSION(klon)                 :: ubwk, vbwk      ! ubwk et vbwk sont les vitesses moyennes dans les poches
REAL, DIMENSION(klon)                 :: weilambda, U10mMOD

INTEGER :: nwb
INTEGER :: i, nmc, kwb
REAL    :: pi, pdfu
REAL    :: auxreal, kref
REAL    :: ray, ray2, theta,rr, xx, yy, zz

pi=2.*acos(0.)
ray=7000.
!ktwk=0.5
!ktth=2.
!kzth=1.
kref=3
nwb=nsurfwind

Do i=1,klon
    ubwk(i) = zu10m(i)
    vbwk(i) = zv10m(i)
ENDDO

DO i=1,klon
     U10mMOD(i)=MAX(woff,sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i)))    
     !U10mMOD(i)=sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
ENDDO

! a enlever 
!call histogram0(0., 20., nbin, hist)
!call histogram0(-20., 20., nbin1, histx)
!call histogram0(-20., 20., nbin1, histy)
! Utilisation du vent moyen dans la maille
IF (iflag_surf_wind == 0) THEN
    !U10mMOD=sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
    IF (nsurfwind /= 1 ) THEN
            STOP 'Si iflag_surf_wind=0, nsurfwind=1'
    ENDIF
    DO i=1,klon
         wind10ms(i,1)=U10mMOD(i)
         probu(i,1)=1.
    ENDDO

!-----------------------------------------------------------------------------
ELSE IF (iflag_surf_wind == 1) THEN ! Weibull
!-----------------------------------------------------------------------------

    DO i=1, klon
        DO nmc=1, nsurfwind
             ! Utilisation de la distribution de weibull
             !U10mMOD=sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
             auxreal=1.+1./kref
             weilambda(i) = U10mMOD(i)/exp(auxreal*log(auxreal)-auxreal &
             - 0.5*log(auxreal/(2.*pi))+1./(12.*auxreal) &
             -1./(360.*(auxreal**3.))+1./(1260.*(auxreal**5.)))
             kwb=nmc
             wind10ms(i,kwb)=kwb*2.*U10mMOD(i)/nwb
    !        pdfu=(kref/U10mMOD(i))*(wind10ms(kwb)/U10mMOD(i))**(kref-1) &
    !        *exp(-(wind10ms(kwb)/U10mMOD(i))**kref)
             pdfu=(kref/weilambda(i))*(wind10ms(i,kwb)/weilambda(i))**(kref-1) &
             *exp(-(wind10ms(i,kwb)/weilambda(i))**kref)
    !        !print *,'JEdbg  U10mMOD(i) weilambda  ',U10mMOD(i),weilambda
             !JE20141205>>
             probu(i,kwb)=pdfu*2.*U10mMOD(i)/nwb
        ENDDO                          
    ENDDO

!-----------------------------------------------------------------------------
ELSE  ! Monte Carlo
!-----------------------------------------------------------------------------

    DO i=1, klon
        DO nmc=1, nsurfwind
            ! Utilisation de la distribution du vent a l interieur et a l exterieur des poches 
            call Random_number(zz)     ! tirage uniforme entre 0 et 1.
            IF (zz <= sigmaw(i)) THEN    ! quand on est a l interieur de la poche
            !IF (zz <= 1.) THEN    ! test pour tourner uniquement avec le modele de vent dans les poches
 
                  call Random_number(xx)   
                  ray2=xx*ray**2
                  call Random_number(yy)
                  theta=yy*2.*pi
                  rr=sqrt(ray2)
                  xp(i,nmc)=rr*cos(theta)
                  yp(i,nmc)=rr*sin(theta)
  

                  ! ENDIFle vent dans la poche = le module du vent moyen dans la poche + un vent radial 
                  vwx(i,nmc) = ubwk(i) + xp(i,nmc)*cstar(i)/ray
                  vwy(i,nmc) = vbwk(i) + yp(i,nmc)*cstar(i)/ray
                  vw(i,nmc) = sqrt(vwx(i,nmc)**2 + vwy(i,nmc)**2)

                  ! On relie la variance au module du vent au carree (sigma ^ 2 = k || v || ^ 2)
                  !sigma_wk(i,nmc) =  ktwk*(vw(i,nmc))
                  sigma_wk(i,nmc) =  surf_wind_ktwake*(vw(i,nmc))

                  ! tirage du vent turbulent vt
                  call randong(xx,yy,pi)
                  vtx(i,nmc) = sigma_wk(i,nmc)*xx
                  vty(i,nmc) = sigma_wk(i,nmc)*yy
                 ! print*,'ZZZZZZZZZZZZZZZZZZZZ xx=',xx
                  ! vent total = vent dans la poche (vw) + le vent turbulent(vt)
                  windx(i,nmc) = vwx(i,nmc) + vtx(i,nmc)
                  windy(i,nmc) = vwy(i,nmc) + vty(i,nmc)
                  wind(i,nmc) = sqrt(windx(i,nmc)**2 + windy(i,nmc)**2)
                  wind10ms(i,nmc) = wind(i,nmc)
                  probu(i,nmc) = 1./nsurfwind

            ELSE
                  ! le vent en dehors des poches est donne par le vent moyen
                  !vwx(klon) = zu10m(klon) 
                  !vwy(klon) = zv10m(klon)
                  !vw(klon) = sqrt(vwx(klon)**2 + vwy(klon)**2)

                  !sigma_th(i,nmc) = sqrt((ktth*ustar(i))**2 + (kzth*wstar(i))**2)  ! a voir
                  ! On remplace wstar par sqrt(2*ale_bl)
                  sigma_th(i,nmc) = sqrt((surf_wind_kttherm*ustar(i))**2 + (surf_wind_kztherm*sqrt(2*ale_bl(i)))**2)                  

                  ! tirage du vent turbulent vt
                  call randong(xx,yy,pi)
                  vtx(i,nmc) = sigma_th(i,nmc)*xx
                  vty(i,nmc) = sigma_th(i,nmc)*yy
  
                  ! vent total en dehors des poches = le vent moyen + le vent turbulent(vt) 
                  windx(i,nmc) = zu10m(i) + vtx(i,nmc)
                  windy(i,nmc) = zv10m(i) + vty(i,nmc)
                  wind(i,nmc) = sqrt(windx(i,nmc)**2 + windy(i,nmc)**2)
                  wind10ms(i,nmc) = wind(i,nmc)
                  probu(i,nmc) = 1./nsurfwind
                  ! print*, 'wind10ms', wind10ms(i,nmc)         
            ENDIF
    ! print*,'WWWWWWWWWWWWWWWWWWWW wind10ms=',wind10ms(i,nmc)
    ! enlver      
    ! call histogram(wind(i,nmc), 0., 20., nbin, hist)
    ! call histogram(windx(i,nmc), -20., 20., nbin1, histx)
    ! call histogram(windy(i,nmc), -20., 20., nbin1, histy)
       ENDDO
   ENDDO
ENDIF
! a enlever
!Do i=1,nbin
!  write(10,*) hist(1,i), hist(2,i)
!enddo

!Do i=1,nbin1
!  write(11,*) histx(1,i), histx(2,i)
!enddo

!Do i=1,nbin1
!  write(12,*) histx(1,i), histy(2,i)
!enddo

RETURN
END SUBROUTINE surf_wind

SUBROUTINE randong(zzz1,zzz2,pi)   ! tirage sur une gaussienne selon BOX-MULLER

implicit none

real      zzz1, zzz2
real      u1, u2
real      pi

! tirage entre 0 et 1 selon la loi uniforme
call RANDOM_NUMBER(u1)
call RANDOM_NUMBER(u2)

! transformation de u1 et u2 en une distribution gaussienne centree reduite sur x et y
zzz1=sqrt(-2*log(u1))*cos(2.*pi*u2)
zzz2=sqrt(-2*log(u1))*sin(2.*pi*u2)

RETURN
END SUBROUTINE randong
END MODULE lmdz_surf_wind
