source: trunk/LMDZ.TITAN/libf/phytitan/sources.F @ 1243

Last change on this file since 1243 was 808, checked in by slebonnois, 12 years ago

SL: Many changes for VENUS (related to newstart) and TITAN (related to clouds). Please read DOC/chantiers/commit_importants.log (cf v808).
File size: 6.7 KB
RevLine 
[175]1         SUBROUTINE sources(ngrid,nlay,
2     $                      ptimestep,pz0,pu,pv,
3     $                      pplev,pzlay,pzlev,
4     $                      gaz1,gaz2,gaz3,
[474]5     $                      ptsrf,evapch4,reserv)
[175]6
7c=======================================================================
8c     Calcul des flux aux interfaces pour les sources
9c     CH4 a la surface
10c     Production de C2H6 en haut du modele.
11c     Production de C2H2 en haut du modele.
12c     
13c     NOTE :
14c     Les gaz ont la tete en haut.
15c     ils ne suivent pas la meme convention que muphys :
16c     (1 -> sol  / klev = haut du modele)
17c=======================================================================
18
19c-----------------------------------------------------------------------
20c   declarations:
21c   -------------
22
23         use dimphy
24         IMPLICIT NONE
25#include "dimensions.h"
26#include "microtab.h"
27#include "YOMCST.h"
28c
29c   arguments:
30c   ----------
31
32         INTEGER ngrid,nlay,nq,ihor
33         REAL ptimestep
34         REAL pplev(ngrid,nlay+1)
35         REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
36         REAL pu(ngrid),pv(ngrid)
37         REAL gaz1(ngrid,nlay),gaz2(ngrid,nlay),gaz3(ngrid,nlay)
38         REAL ptsrf(ngrid)
[474]39         REAL evapch4(ngrid)
[175]40c
41c   local:
42c   ------
43 
44         INTEGER ilev,ig,ilay,nlev,k,inch4,inc2h6
45
46         REAL zgz1(klon,klev),zgz2(klon,klev),zgz3(klon,klev)
47         REAL zcdv(klon),zu2,pz0
48         REAL xmair,gg,zrho,ws,ch,qch4,flux
[808]49         REAL effg               ! effg est une fonction(z), z en m.
[175]50         REAL xmuair
51         REAL zmem,zmem2,zmem3
52         REAL prodc2h6,prodc2h2
53         real reserv(ngrid),restemp,drestemp
[474]54         REAL zevapch4
[175]55
56         real umin
57         data umin/1.e-12/
58         save umin
59c
[474]60c
[175]61c-----------------------------------------------------------------------
62c   initialisations:
63c   -----------------
64
65         nlev=nlay+1
66
67         if(nlay.ne.klev) THEN
68           PRINT*,'STOP dans sources.F'
69           PRINT*,'probleme de dimensions :'
70           PRINT*,'nlay  =',nlay
71           PRINT*,'klev  =',klev
72           STOP
73         endif
74
75         IF(ngrid.NE.klon) THEN
76           PRINT*,'STOP dans sources.F'
77           PRINT*,'probleme de dimensions :'
78           PRINT*,'ngrid  =',ngrid
79           PRINT*,'klon  =',klon
80           STOP
81         ENDIF
82
83         zgz1 = gaz1
84         zgz2 = gaz2
85         zgz3 = gaz3
86
[474]87         evapch4 = 0.
88
[175]89c-----------------------------------------------------------------------
90c     2. calcul de  cd :
91c     ----------------
92c
93         DO ig=1,ngrid
94           zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin
[808]95           zcdv(ig)=pz0*(sqrt(zu2))
[175]96c           write(99,'(I4,3(ES24.17,1X))') ig,
[808]97c     &     pz0,zu2,(sqrt(zu2))
[175]98         ENDDO
99c          write(99,*) ""
100
101c-----------------------------------------------------------------------
102c    4. Conditions aux limites pour CH4 et C2H6
103c    -------------------------------------------
104c
105
106*   Conditions CH4
107         DO ig=1,ngrid
[474]108           zevapch4=0.
[175]109           restemp=0.
[808]110           gg=effg(pzlay(ig,1))
[175]111           zrho=(pplev(ig,1)-pplev(ig,2))/gg
112           zrho=zrho/(pzlev(ig,2)-pzlev(ig,1))
113           ws=sqrt(pu(ig)**2.+pv(ig)**2.)*(10./pzlay(ig,1))**0.2
114           ch=1.5*sqrt(zcdv(ig))
115           call ch4sat(ptsrf(ig),pplev(ig,1),qch4) ! qch4=kg/kg
116           qch4=qch4*0.50 ! ici on impose 50% d'humidité au sol
117       
118           if(reserv(ig).le. 1.e-10 ) then 
119             flux=0.
120             reserv(ig)=1.e-10
121           else
122             flux=zrho*ch*ws
123             flux=flux*0.1 ! fraction occupée par les lacs
124           endif
125
126           zmem=zgz1(ig,1)
127           zgz1(ig,1)=(zgz1(ig,1)+flux*ptimestep*qch4*28./16.)
128     &                /(1.+flux*ptimestep)
129
[808]130           gg=effg(pzlay(ig,1))
[175]131           xmair=(pplev(ig,1)-pplev(ig,1+1))/gg
132           xmair=xmair/(pzlev(ig,1+1)-pzlev(ig,1))
133           xmuair=28.!*(1.-zmem)+zmem*16.
134
135           drestemp = - (zgz1(ig,1)-zmem)*xmair ! en m3/m2=m
136     &     *(pzlev(ig,2)-pzlev(ig,1))*16./xmuair/425.
137
138c           ici on peut fixer un seuil sur drestemp
139c           (ie limiter l'echange atm/surface)
140     
141           restemp=reserv(ig) +drestemp
142
143           IF (restemp.ge.0.) THEN
144             reserv(ig) = reserv(ig) + drestemp
[474]145             zevapch4   = zevapch4   + drestemp
[175]146           ELSE
147*          Il n'y a pas suffisamment de méthane; on re-évalue le flux d'évaporation
148*          Quelle nouvelle concentration zgz1(ig,1) atteint-on en évaporant tout ?
[474]149             zgz1(ig,1)= reserv(ig)/(xmair*(pzlev(ig,1+1)-pzlev(ig,1))
[175]150     &                 *16./xmuair/425.)+zmem
[474]151             zevapch4  = zevapch4-reserv(ig)
[175]152
153             if(reserv(ig).eq.0.)
154     &       print*,'assechement du sol en ig=', ig,reserv(ig),flux
155
156             reserv(ig)=0.  ! on a tout évaporé...
157           ENDIF
158c         
[808]159           evapch4(ig) = zevapch4  ! < 0 si volume évaporé (m3/m2)
[175]160
161         ENDDO
[474]162
[175]163*   Conditions C2H6
164
165         prodc2h6=6.e-12/5. ! kg/m2/s
166 
167         IF (1.EQ.1) THEN
168           DO ig=1,ngrid
169             DO ilev=nlay,nlay-4,-1
170*            calcule de zrho (kg/m3) pour la couche...
[808]171               gg=effg(pzlay(ig,ilev))
[175]172               zrho=(pplev(ig,ilev)-pplev(ig,ilev+1))/gg
173               zrho=zrho/(pzlev(ig,ilev+1)-pzlev(ig,ilev))
174
175*              passage taux de production --> Dx_c2h6 a rajouter au niveau ilev
176               zmem2=zgz2(ig,ilev)
177               zgz2(ig,ilev)=zgz2(ig,ilev)+
178     &         prodc2h6*ptimestep/
179     &         abs(pzlev(ig,ilev+1)-pzlev(ig,ilev))    !  kg/m3/s
180     &         /zrho*28./30.
181
182             ENDDO
183           ENDDO
184
185         ELSE
186
187           DO ig=1,ngrid
188             DO ilev=nlay,nlay-8,-1
189               zgz2(ig,ilev)=1.2e-5
190             ENDDO
191           ENDDO
192
193         ENDIF  ! (fin 1.eq.0)
194
195*-------------------------------------
196*   Conditions C2H2
197
198         prodc2h2=1.6e-12/5. ! kg/m2/s
199
200         IF(1  .EQ.  1) THEN
201
202           DO ig=1,ngrid
203             DO ilev=nlay,nlay-4,-1
204*            calcule de zrho (kg/m3) pour la couche...
[808]205                gg=effg(pzlay(ig,ilev))
[175]206                zrho=(pplev(ig,ilev)-pplev(ig,ilev+1))/gg
207                zrho=zrho/(pzlev(ig,ilev+1)-pzlev(ig,ilev))
208
209*            passage taux de production --> Dx_c2h2 a rajouter au niveau ilev
210                zmem3=zgz3(ig,ilev)
211                zgz3(ig,ilev)=zgz3(ig,ilev)+
212     &          (prodc2h2)*ptimestep/
213     &          abs(pzlev(ig,ilev+1)-pzlev(ig,ilev))  !  kg/m3/s
214     &          /zrho*28./26.
215
216             ENDDO
217
218           ENDDO
219
220         ENDIF
221
222c-----------------------------------------------------------------------
223         DO ig=1,ngrid
224           DO ilev=1,nlay
225             gaz1(ig,ilev)=zgz1(ig,ilev)
226             gaz2(ig,ilev)=zgz2(ig,ilev)
227             gaz3(ig,ilev)=zgz3(ig,ilev)
228           ENDDO
229         ENDDO
230
231        RETURN
232        END
Note: See TracBrowser for help on using the repository browser.