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

Last change on this file since 3094 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

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