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