Index: /LMDZ6/trunk/libf/phylmd/lmdz_surf_wind.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_surf_wind.f90	(revision 5445)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_surf_wind.f90	(revision 5445)
@@ -0,0 +1,189 @@
+MODULE lmdz_surf_wind
+        CONTAINS 
+
+SUBROUTINE surf_wind(klon,nsrfwnd,zu10m,zv10m,sigmaw,cstar,ustar,wstar,wind10ms,probu)
+
+USE lmdz_surf_wind_ini, ONLY : iflag_surf_wind
+
+IMPLICIT NONE
+INTEGER, INTENT(IN)                :: nsrfwnd, 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, wstar
+REAL, DIMENSION(klon,nsrfwnd), INTENT(OUT)         :: wind10ms, probu
+
+
+REAL, DIMENSION(klon,nsrfwnd)         :: sigma_th, sigma_wk
+REAL, DIMENSION(klon,nsrfwnd)         :: xp, yp, zz
+REAL, DIMENSION(klon,nsrfwnd)         :: vwx, vwy, vw
+REAL, DIMENSION(klon,nsrfwnd)         :: vtx, vty
+REAL, DIMENSION(klon,nsrfwnd)         :: 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
+REAL    :: ktwk, ktth, kzth
+
+print*,'LLLLLLLLLLLLLLLLLLLLL nsrfwnd=',nsrfwnd
+pi=2.*acos(0.)
+ray=7000.
+ktwk=0.5
+ktth=2.
+kzth=1.
+kref=3
+nwb=nsrfwnd
+
+ubwk(klon) = zu10m(klon)
+vbwk(klon) = zv10m(klon)
+
+DO i=1,klon
+    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 (nsrfwnd /= 1 ) THEN
+            STOP 'Si iflag_surf_wind=0, nsrfwnd=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, nsrfwnd
+             ! 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, nsrfwnd
+            ! 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 (ALL(zz <= sigmaw(klon))) THEN    ! quand on est a l interieur de la poche
+ 
+                  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))
+
+                  ! 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
+
+                  ! 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) = wind(i,nmc)/nsrfwnd
+
+            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
+                  sigma_th(i,nmc) = 1.8
+
+                  ! 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) = wind(i,nmc)/nsrfwnd
+                  ! print*, 'wind10ms', wind10ms(i,nmc)         
+            ENDIF
+    ! 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
Index: /LMDZ6/trunk/libf/phylmd/lmdz_surf_wind_ini.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_surf_wind_ini.f90	(revision 5445)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_surf_wind_ini.f90	(revision 5445)
@@ -0,0 +1,48 @@
+MODULE lmdz_surf_wind_ini
+IMPLICIT NONE 
+save
+
+  ! ============================================================================
+  ! Declaration de variables
+  ! -------------------------------------------------------------------------
+
+   integer, protected :: lunout
+   integer, protected :: iflag_surf_wind=0
+   !$OMP THREADPRIVATE(lunout,iflag_surf_wind,nsrfwnd)
+
+CONTAINS
+
+  ! =========================================================================
+  SUBROUTINE surf_wind_ini(klon,lunout_i)
+  ! =========================================================================
+
+  ! **************************************************************
+  ! *
+  ! WAKE                                                        *
+  ! retour a un Pupper fixe                                *
+  ! *
+  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
+  ! modified by :   ROEHRIG Romain        01/29/2007            *
+  ! **************************************************************
+
+  ! -------------------------------------------------------------------------
+  ! Initialisations
+  ! -------------------------------------------------------------------------
+
+
+  USE ioipsl_getin_p_mod, ONLY : getin_p
+  integer, intent(in) :: klon,lunout_i
+
+  lunout=lunout_i
+  CALL getin_p('iflag_surf_wind',iflag_surf_wind)
+
+  write(lunout,*) 'Initialisation wind10m'
+  write(lunout,*) 'lmdz_surf_wind_ini, iflag_surf_wind=',iflag_surf_wind
+
+ RETURN
+
+END SUBROUTINE surf_wind_ini
+
+
+
+END MODULE lmdz_surf_wind_ini
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5444)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5445)
@@ -80,4 +80,6 @@
     USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
     USE lmdz_wake_ini, ONLY : wake_ini
+    USE lmdz_surf_wind_ini, ONLY : surf_wind_ini, iflag_surf_wind
+    USE lmdz_surf_wind, ONLY : surf_wind
     USE yamada_ini_mod, ONLY : yamada_ini
     USE lmdz_atke_turbulence_ini, ONLY : atke_ini
@@ -1263,4 +1265,13 @@
     CHARACTER(len=512) :: namelist_ecrad_file
 
+
+    ! Subgrid scale wind :
+    ! Need to be allocatable/save because the number of bin is not known (provided by surf_wind_ini)
+    integer, save :: nsrfwnd=1
+    real, dimension(:,:), allocatable, save :: surf_wind_value, surf_wind_proba ! module and probability of sugrdi wind wind sample
+    !$OMP THREADPRIVATE(nsrfwnd,surf_wind_value, surf_wind_proba)
+    
+
+
     !======================================================================!
     ! Bifurcation vers un nouveau moniteur physique pour experimenter      !
@@ -1818,4 +1829,9 @@
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       CALL surf_wind_ini(klon,lunout)
+       CALL getin_p('nsrfwnd',nsrfwnd)
+       allocate(surf_wind_value(klon,nsrfwnd),surf_wind_proba(klon,nsrfwnd))
+    
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        CALL wake_ini(rg,rd,rv,prt_level)
@@ -3802,4 +3818,8 @@
     !
     !===================================================================
+    ! Computation of subrgid scale near-surface wind distribution
+    call surf_wind(klon,nsrfwnd,u10m,v10m,wake_s,wake_Cstar,ustar,wstar,surf_wind_value,surf_wind_proba)
+
+    !===================================================================
     ! Computation of ratqs, the width (normalized) of the subrid scale
     ! water distribution
