source: LMDZ6/trunk/libf/phylmd/lmdz_surf_wind.f90 @ 5715

Last change on this file since 5715 was 5635, checked in by fhourdin, 8 weeks ago

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