source: LMDZ6/branches/contrails/libf/phylmd/lmdz_surf_wind.f90 @ 5450

Last change on this file since 5450 was 5450, checked in by aborella, 7 hours ago

Merge with trunk
File size: 6.5 KB
Line 
1MODULE lmdz_surf_wind
2        CONTAINS
3
4SUBROUTINE surf_wind(klon,nsrfwnd,zu10m,zv10m,sigmaw,cstar,ustar,wstar,wind10ms,probu)
5
6USE lmdz_surf_wind_ini, ONLY : iflag_surf_wind
7
8IMPLICIT NONE
9INTEGER, INTENT(IN)                :: nsrfwnd, klon
10REAL, DIMENSION(klon), INTENT(IN)  :: zu10m, zv10m
11REAL, DIMENSION(klon), INTENT(IN)  :: cstar
12REAL, DIMENSION(klon), INTENT(IN)  :: sigmaw
13REAL, DIMENSION(klon), INTENT(IN)  :: ustar, wstar
14REAL, DIMENSION(klon,nsrfwnd), INTENT(OUT)         :: wind10ms, probu
15
16
17REAL, DIMENSION(klon,nsrfwnd)         :: sigma_th, sigma_wk
18REAL, DIMENSION(klon,nsrfwnd)         :: xp, yp, zz
19REAL, DIMENSION(klon,nsrfwnd)         :: vwx, vwy, vw
20REAL, DIMENSION(klon,nsrfwnd)         :: vtx, vty
21REAL, DIMENSION(klon,nsrfwnd)         :: windx, windy, wind
22REAL, DIMENSION(klon)                 :: ubwk, vbwk      ! ubwk et vbwk sont les vitesses moyennes dans les poches
23REAL, DIMENSION(klon)                 :: weilambda, U10mMOD
24
25INTEGER :: nwb
26INTEGER :: i, nmc, kwb
27REAL    :: pi, pdfu
28REAL    :: auxreal, kref
29REAL    :: ray, ray2, theta,rr, xx, yy
30REAL    :: ktwk, ktth, kzth
31
32print*,'LLLLLLLLLLLLLLLLLLLLL nsrfwnd=',nsrfwnd
33pi=2.*acos(0.)
34ray=7000.
35ktwk=0.5
36ktth=2.
37kzth=1.
38kref=3
39nwb=nsrfwnd
40
41ubwk(klon) = zu10m(klon)
42vbwk(klon) = zv10m(klon)
43
44DO i=1,klon
45    U10mMOD(i)=sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
46ENDDO
47
48! a enlever
49!call histogram0(0., 20., nbin, hist)
50!call histogram0(-20., 20., nbin1, histx)
51!call histogram0(-20., 20., nbin1, histy)
52! Utilisation du vent moyen dans la maille
53IF (iflag_surf_wind == 0) THEN
54    !U10mMOD=sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
55    IF (nsrfwnd /= 1 ) THEN
56            STOP 'Si iflag_surf_wind=0, nsrfwnd=1'
57    ENDIF
58    DO i=1,klon
59         wind10ms(i,1)=U10mMOD(i)
60         probu(i,1)=1.
61    ENDDO
62
63!-----------------------------------------------------------------------------
64ELSE IF (iflag_surf_wind == 1) THEN ! Weibull
65!-----------------------------------------------------------------------------
66
67    DO i=1, klon
68        DO nmc=1, nsrfwnd
69             ! Utilisation de la distribution de weibull
70             !U10mMOD=sqrt(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
71             auxreal=1.+1./kref
72             weilambda(i) = U10mMOD(i)/exp(auxreal*log(auxreal)-auxreal &
73             - 0.5*log(auxreal/(2.*pi))+1./(12.*auxreal) &
74             -1./(360.*(auxreal**3.))+1./(1260.*(auxreal**5.)))
75             kwb=nmc
76             wind10ms(i,kwb)=kwb*2.*U10mMOD(i)/nwb
77    !        pdfu=(kref/U10mMOD(i))*(wind10ms(kwb)/U10mMOD(i))**(kref-1) &
78    !        *exp(-(wind10ms(kwb)/U10mMOD(i))**kref)
79             pdfu=(kref/weilambda(i))*(wind10ms(i,kwb)/weilambda(i))**(kref-1) &
80             *exp(-(wind10ms(i,kwb)/weilambda(i))**kref)
81    !        !print *,'JEdbg  U10mMOD(i) weilambda  ',U10mMOD(i),weilambda
82             !JE20141205>>
83             probu(i,kwb)=pdfu*2.*U10mMOD(i)/nwb
84        ENDDO                         
85    ENDDO
86
87!-----------------------------------------------------------------------------
88ELSE  ! Monte Carlo
89!-----------------------------------------------------------------------------
90
91    DO i=1, klon
92        DO nmc=1, nsrfwnd
93            ! Utilisation de la distribution du vent a l interieur et a l exterieur des poches
94            call Random_number(zz)     ! tirage uniforme entre 0 et 1.
95            IF (ALL(zz <= sigmaw(klon))) THEN    ! quand on est a l interieur de la poche
96 
97                  call Random_number(xx)   
98                  ray2=xx*ray**2
99                  call Random_number(yy)
100                  theta=yy*2.*pi
101                  rr=sqrt(ray2)
102                  xp(i,nmc)=rr*cos(theta)
103                  yp(i,nmc)=rr*sin(theta)
104 
105
106                  ! ENDIFle vent dans la poche = le module du vent moyen dans la poche + un vent radial
107                  vwx(i,nmc) = ubwk(i) + xp(i,nmc)*cstar(i)/ray
108                  vwy(i,nmc) = vbwk(i) + yp(i,nmc)*cstar(i)/ray
109                  vw(i,nmc) = sqrt(vwx(i,nmc)**2 + vwy(i,nmc)**2)
110
111                  ! On relie la variance au module du vent au carree (sigma ^ 2 = k || v || ^ 2)
112                  sigma_wk(i,nmc) =  ktwk*(vw(i,nmc))
113
114                  ! tirage du vent turbulent vt
115                  call randong(xx,yy,pi)
116                  vtx(i,nmc) = sigma_wk(i,nmc)*xx
117                  vty(i,nmc) = sigma_wk(i,nmc)*yy
118
119                  ! vent total = vent dans la poche (vw) + le vent turbulent(vt)
120                  windx(i,nmc) = vwx(i,nmc) + vtx(i,nmc)
121                  windy(i,nmc) = vwy(i,nmc) + vty(i,nmc)
122                  wind(i,nmc) = sqrt(windx(i,nmc)**2 + windy(i,nmc)**2)
123                  wind10ms(i,nmc) = wind(i,nmc)
124                  probu(i,nmc) = wind(i,nmc)/nsrfwnd
125
126            ELSE
127                  ! le vent en dehors des poches est donne par le vent moyen
128                  !vwx(klon) = zu10m(klon)
129                  !vwy(klon) = zv10m(klon)
130                  !vw(klon) = sqrt(vwx(klon)**2 + vwy(klon)**2)
131
132                  !sigma_th(i,nmc) = sqrt((ktth*ustar(i))**2 + (kzth*wstar(i))**2)  ! a voir
133                  sigma_th(i,nmc) = 1.8
134
135                  ! tirage du vent turbulent vt
136                  call randong(xx,yy,pi)
137                  vtx(i,nmc) = sigma_th(i,nmc)*xx
138                  vty(i,nmc) = sigma_th(i,nmc)*yy
139 
140                  ! vent total en dehors des poches = le vent moyen + le vent turbulent(vt)
141                  windx(i,nmc) = zu10m(i) + vtx(i,nmc)
142                  windy(i,nmc) = zv10m(i) + vty(i,nmc)
143                  wind(i,nmc) = sqrt(windx(i,nmc)**2 + windy(i,nmc)**2)
144                  wind10ms(i,nmc) = wind(i,nmc)
145                  probu(i,nmc) = wind(i,nmc)/nsrfwnd
146                  ! print*, 'wind10ms', wind10ms(i,nmc)         
147            ENDIF
148    ! enlver     
149    !call histogram(wind(i,nmc), 0., 20., nbin, hist)
150    !call histogram(windx(i,nmc), -20., 20., nbin1, histx)
151    !call histogram(windy(i,nmc), -20., 20., nbin1, histy)
152       ENDDO
153   ENDDO
154ENDIF
155! a enlever
156!Do i=1,nbin
157!  write(10,*) hist(1,i), hist(2,i)
158!enddo
159
160!Do i=1,nbin1
161!  write(11,*) histx(1,i), histx(2,i)
162!enddo
163
164!Do i=1,nbin1
165!  write(12,*) histx(1,i), histy(2,i)
166!enddo
167
168RETURN
169END SUBROUTINE surf_wind
170
171SUBROUTINE randong(zzz1,zzz2,pi)   ! tirage sur une gaussienne selon BOX-MULLER
172
173implicit none
174
175real      zzz1, zzz2
176real      u1, u2
177real      pi
178
179! tirage entre 0 et 1 selon la loi uniforme
180call RANDOM_NUMBER(u1)
181call RANDOM_NUMBER(u2)
182
183! transformation de u1 et u2 en une distribution gaussienne centree reduite sur x et y
184zzz1=sqrt(-2*log(u1))*cos(2.*pi*u2)
185zzz2=sqrt(-2*log(u1))*sin(2.*pi*u2)
186
187RETURN
188END SUBROUTINE randong
189END MODULE lmdz_surf_wind
Note: See TracBrowser for help on using the repository browser.