source: trunk/LMDZ.TITAN/libf/phytitan/cnuages3D.F @ 1242

Last change on this file since 1242 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: 13.6 KB
RevLine 
[175]1         subroutine cnuages(
2     &   tq,tqc1,tqc2,tqc3,tqcn,gaz1,gaz2,gaz3,   ! aerosol/glace/gas
3     &   ddt)
4
5c
6c
7c  °
8c              °
9c  SERT A APPELE LA ROUTINE MICROPHYSIQUE DES NUAGES
10c
11c     °
12c  ICI ON NE FAIT QUE LA NUCLEATION/CONDENSATION
13c  ET GESTION°DES NOYAUX. LA SEDIMENTATION EST DANS
14c  SNUAGES.F
15c            °   °
16c
17c        °
18c                °
19c              °   °
20c                     
21c                  °   °
22c                   ° 
23c                    ° °
24c                   ° \|/  °
25c                    (@ @)°
26c-----------------oOo--O--oOo--------------------------
27c                       
28c
29c
30c Interface entre physiq.F et les routines n_<nom_compose>.F
31c
32c Date: 3 Nov 2003
33c
34c   EN ENTREE/SORTIE DE LA ROUTINE
35c ------------------------------------
36c
37c     Les aerosols, noyaux (tq,tqcn) sont en nbre/m^2 dans la colonne
38c     Les  condensats (tqc1,tqc2)   sont en volume/m^2 dans la colonne
39c     Le gaz (gaz1, gaz2) est en fraction molaire
40c
41c   EN APPEL DES ROUTINES NUAGES
42c ------------------------------------
43c
44c     Les aerosols et noyaux doivent etre en nombre /kg d'air
45c     Les condensats doivent etre en volume / kg d'air
46c     Le gaz en kg/kg d'air
47c
48c   LES TENDANCES ET DIFFERENCES SONT HOMOGENES AUX QUANTITES
49c ------------------------------------------------------------
50c
51c
52c------------------------------------------------------
53
54         use dimphy
55         IMPLICIT NONE
56#include "dimensions.h"
57#include "microtab.h"
58#include "varmuphy.h"
59
60         integer NG1,NG,NL
61         parameter (NG1=1,NG=NG1,NL=llm)
62 
63c*************************************
64c  declaration des variables internes *
65c*************************************
66
67c        INTERNE!        *
68c-----------------------*
69
70         real  tqc1(NG,NL,nrad)
71         real  tqc2(NG,NL,nrad)
72         real  tqc3(NG,NL,nrad)
73         real  tqcn(NG,NL,nrad)
74         real  tq(NG,NL,nrad)
75*
76         real  tdqc1(NG,NL,nrad)
77         real  tdqc2(NG,NL,nrad)
78         real  tdqc3(NG,NL,nrad)
79         real  tdq(NG,NL,nrad,ntype-2+1)
80         real  tdqcn(NG,NL,nrad,ntype-2+1)
81*
82         real  gaz1(NG,NL)
83         real  gaz2(NG,NL)
84         real  gaz3(NG,NL)
85         real  dgaz1(NG,NL)
86         real  dgaz2(NG,NL)
87         real  dgaz3(NG,NL)
88*
89         real  ppch4(NG,NL)
90         real  ppc2h6(NG,NL)
91         real  ppn2(NG,NL)
92*
93         real  pmixch4(NL)
94         real  pmixc2h6(NL)
95         real  pmixn2(NL)
96c   composition initiale estimée (interne)
97         real  xprime1(NG,NL)
98         real  xprime2(NG,NL)
99         real  xprime3(NG,NL)
100c   composition calculée (output)
101         real  x1(NL)
102         real  x2(NL)
103         real  x3(NL)
104c   moyenne "glissante" pondéréee (output + mémoire)
105         real  x1o(NL)
106         real  x2o(NL)
107         real  x3o(NL)
108         real  icefrac(NL)
109         real  dmn2(NL+1)
110
111         real  ppch4t,ppc2h6t,ppn2t
112         real  psatch4,psatc2h6,psatn2
113         real  xprime(3),x(3),frac
114         real  melange
115         real  sum,sum0
116
117*                    RAPPEL: NG=1
118
119         real ddt
120         real masspaer
121
122         common/mixing/x1,x2,x3,icefrac,
123     &         pmixch4,pmixc2h6,pmixn2,
124     &         x1o,x2o,x3o
125
126
127
128c  FORMAT MICRO DES NUAGES
129c------------------------*
130
131         real  especes(NG,NL,3*nrad+1)         
132         real  condens(NG,NL,nrad)         
[808]133         real  gg(NL),xmair
134         real  effg     ! effg est une fonction(z), z en m.
[175]135
136         integer jsup,jinf,h,i,j,k,ndim
137         integer ival1,ival2,ival3
138
139         integer iprem
140
141         save iprem,xprime
142         data iprem/0/
143
144         ndim=3*nrad+1
145
[808]146         do j=1,NL
147           gg(j)=effg(z(j))
148         enddo
[175]149
150*********************************************
151*********************************************
152*  Appel de la condensation du methane
153*********************************************
154*********************************************
155
156
157
158*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
159*                   Bilan avant sur le methane                         *
160*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
161
162         do i=1,ng1  ! ng1=1 !!
163           do j=1,NL
164
165*  RAZ des delta d'especes communes.
166*----------------------------------
167             do k=1,nrad
168               tdqcn(i,j,k,1)= 0.
169               tdq( i,j,k,1) = 0.
170             enddo
171
[808]172             xmair=(pb(j+1)-pb(j))/gg(j)/dzb(j)
[175]173
174             do k=1,nrad
175               especes(i,j,k)=tq(i,j,k)         /xmair       ! aerosols, noyaux,
176               especes(i,j,k+nrad)=tqc1(i,j,k)  /xmair       ! methane condense,
177               condens(i,j,k)=(tqc2(i,j,k)+tqc3(i,j,k))/xmair! autre(s) condensat(s)
178               especes(i,j,k+2*nrad)=tqcn(i,j,k)/xmair       ! nombre/kg &  volume/kg
179             enddo
180             especes(i,j,3*nrad+1)=gaz1(i,j)*mch4/mair       ! methane gazeux kg/kg
181           enddo
182         enddo
183
184
1851001   format(7(1x,e12.6),' avN2CH4C2H6')
1861003   format(7(1x,e12.6),' miN2CH4C2H6')
1871002   format(7(1x,e12.6),' apN2CH4C2H6')
188*
189         call n_methane(ng1,ndim,nrad,ddt,
190     &                  p,t,r_e,especes,condens)
191
192
193*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
194*                   Bilan apres sur le methane                         *
195*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
196
197         do i=1,ng1
198           do j=1,NL
199
[808]200             xmair=(pb(j+1)-pb(j))/gg(j)/dzb(j)
[175]201
202*  ici ce sont les tendances a sortir de nuages.F pour le methane....
203*-------------------------------------------------------------------
204             sum=0.   
205             do k=1,nrad
206               tdqc1(i,j,k)=(especes(i,j,k+nrad)*xmair-tqc1(i,j,k) )
207               sum =sum+tdqc1(i,j,k)/xmair*rhoi_ch4
208               tqc1(i,j,k) = especes(i,j,k+nrad)*xmair
209             enddo
210
211             dgaz1(i,j)= especes(i,j,3*nrad+1)*mair/mch4-gaz1(i,j)
212             gaz1(i,j)= especes(i,j,3*nrad+1)*mair/mch4
213
214c         dgaz1(i,j)=-sum*xmuair/16.
215c          gaz1(i,j)=gaz1(i,j)+dgaz1(i,j)
216
217*  Premiere tendance sur les variables communes (aerosols et noyaux)
218*------------------------------------------------------------------
219             do k=1,nrad
220               tdqcn(i,j,k,1)=(especes(i,j,k+2*nrad)*xmair-tqcn(i,j,k))
221               tdq( i,j,k,1) =(especes(i,j,k)*xmair       -tq(i,j,k)) 
222             enddo
223           enddo
224         enddo
225
226
227
228*     attention, si il y a  de l'ethane sur les noyaux... il est impossible
229*     de les restituer - en revanche on peut en creer de nouveaux ! !!
230*     Le corrolaire de la condition ci dessus est que si il est impossible de
231*     restituer des noyaux, le nombre d'aeorsols ne peux pas augmenter, il
232*     peut en revanche diminuer
233
234
235*********************************************
236*********************************************
237*  Appel de la condensation de l'ethane
238*********************************************
239*********************************************
240
241         do i=1,ng1
242           do j=1,NL
243             do k=1,nrad
244               tdqcn(i,j,k,2)= 0.
245               tdq( i,j,k,2) = 0.
246             enddo
247
[808]248             xmair=(pb(j+1)-pb(j))/gg(j)/dzb(j)
[175]249
250             do k=1,nrad
251               especes(i,j,k)=tq(i,j,k)/xmair
252               especes(i,j,k+nrad)=tqc2(i,j,k)/xmair           ! ethane condense
253               condens(i,j,k)=(tqc1(i,j,k)+tqc3(i,j,k))/xmair  ! autres condensats
254               especes(i,j,k+2*nrad)=tqcn(i,j,k)/xmair
255             enddo
256
257             especes(i,j,3*nrad+1)=gaz2(i,j)*mc2h6/mair       ! ethane gazeux
258
259           enddo
260         enddo
261
262
263
264         call n_ethane(ng1,ndim,nrad,ddt,
265     &                 p,t,r_e,especes,condens)
266
267         do i=1,ng1
268           do j=1,NL
[808]269             xmair=(pb(j+1)-pb(j))/gg(j)/dzb(j)
[175]270
271*  ici ce sont les tendances a sortir de nuages.F pour l'ethane....
272*-----------------------------------------------------------------
273
274             sum=0.
275             do k=1,nrad
276               tdqc2(i,j,k)=(especes(i,j,k+nrad)*xmair-tqc2(i,j,k) )
277               sum =sum+tdqc2(i,j,k)/xmair*rhoi_c2h6
278               tqc2(i,j,k) = especes(i,j,k+nrad)*xmair
279             enddo
280
281             dgaz2(i,j)=especes(i,j,3*nrad+1)*mair/mc2h6 - gaz2(i,j)
282             gaz2(i,j)=especes(i,j,3*nrad+1)*mair/mc2h6
283
284c         dgaz2(i,j)=-sum*xmuair/30.
285c          gaz2(i,j)=gaz2(i,j)+dgaz2(i,j)
286
287
288*  Deuxieme tendance sur les variables communes (aerosols et noyaux)
289*------------------------------------------------------------------
290   
291             do k=1,nrad
292               tdqcn(i,j,k,2)=(especes(i,j,k+2*nrad)*xmair-tqcn(i,j,k))
293               tdq(i,j,k,2)  =(especes(i,j,k)*xmair       -tq(i,j,k))
294             enddo
295
296           enddo
297         enddo
298
299
300*********************************************
301*********************************************
302*  Appel de la condensation de l'acethylene
303*********************************************
304*********************************************
305
306         do i=1,ng1
307           do j=1,NL
308             do k=1,nrad
309               tdqcn(i,j,k,3)= 0.
310               tdq( i,j,k,3) = 0.
311             enddo
[808]312             xmair=(pb(j+1)-pb(j))/gg(j)/dzb(j)
[175]313
314             do k=1,nrad
315               especes(i,j,k)=tq(i,j,k)/xmair
316               especes(i,j,k+nrad)=tqc3(i,j,k)/xmair           ! acethylene condense
317               condens(i,j,k)=(tqc1(i,j,k)+tqc2(i,j,k))/xmair  ! autres condensats
318               especes(i,j,k+2*nrad)=tqcn(i,j,k)/xmair
319             enddo
320
321             especes(i,j,3*nrad+1)=gaz3(i,j)*mc2h2/mair      ! acethylene gazeux
322
323           enddo
324         enddo
325
326         call n_acethylene(ng1,ndim,nrad,ddt,
327     &                     p,t,r_e,especes,condens)
328
329
330         do i=1,ng1
331           do j=1,NL
332
[808]333             xmair=(pb(j+1)-pb(j))/gg(j)/dzb(j)
[175]334
335*  ici ce sont les tendances a sortir de nuages.F pour l'ethane....
336*-----------------------------------------------------------------
337
338             sum=0.
339             do k=1,nrad
340               tdqc3(i,j,k)=(especes(i,j,k+nrad)*xmair-tqc3(i,j,k) )
341               sum =sum+tdqc3(i,j,k)/xmair*rhoi_c2h2
342               tqc3(i,j,k) = especes(i,j,k+nrad)*xmair
343             enddo
344
345             dgaz3(i,j)=especes(i,j,3*nrad+1)*mair/mc2h2 - gaz3(i,j)
346             gaz3(i,j)=especes(i,j,3*nrad+1)*mair/mc2h2
347
348c         dgaz3(i,j)=-sum*xmuair/26.
349c         gaz3(i,j)=gaz3(i,j)+dgaz3(i,j)
350
351*  Troisieme tendance sur les variables communes (aerosols et noyaux)
352*------------------------------------------------------------------
353
354
355             do k=1,nrad
356               tdqcn(i,j,k,3)=(especes(i,j,k+2*nrad)*xmair-tqcn(i,j,k))
357               tdq(i,j,k,3)  =(especes(i,j,k)*xmair       -tq(i,j,k))
358             enddo
359
360           enddo
361         enddo
362
363
364*  FIN DES APPELS DE NUAGES ET BILAN DES TENDANCES...
365*------------------------------------------------------------------
366
367
368         do i=1,ng1
369           do j=1,NL
370             do k=1,nrad
371
372*          Ici on test l'activité nuageuse : si on a l'association
373*          tdqcX(i,j,k) = 0 et tqcX(i,j,k) = 0 alors ivalX reste à 0 (pas d'ativité)
374*                                              sinon ivalX passe à 1 (activité)
375*------------------------------------------------------------------------------------
376
377               ival1=0
378               ival2=0
379               ival3=0
380
381               if(tdqc1(i,j,k).ne.0. .or. tqc1(i,j,k).gt.0.)  ival1=1
382               if(tdqc2(i,j,k).ne.0. .or. tqc2(i,j,k).gt.0.)  ival2=1
383               if(tdqc3(i,j,k).ne.0. .or. tqc3(i,j,k).gt.0.)  ival3=1
384
385*          Ici on definit la tendances des noyaux en faisant deux choses:
386*          -1 On ecarte les cas  tdqcn(i,j,k,X)=0 si ils sont associés à une
387*          absence d'activité nuageuse de l'espèce (tdqcX(i,j,k)=0.)
388*          -2 Sélectionne la tendance la plus élevée. Si aucune activité nuageuse
389*           n'exits dans cette case (ivalX=0 pour les 3 especes), alors on
390*           retrouve la valeur -1.e40 que l'on mets alors à 0.
391*----------------------------------------------------------------------
392c23456789012345678901234567890123456789012345678901234567890123456789012
393
394               tdqcn(i,j,k,ntype-1)=-1.e40                      ! plus petite valeur possible
395
396               if(ival1.eq.1.and.
397     &         tdqcn(i,j,k,1).ge.tdqcn(i,j,k,ntype-1)) ! Si activité de l'espece 1
398     &         tdqcn(i,j,k,ntype-1)=tdqcn(i,j,k,1)
399               if(ival2.eq.1.and.
400     &         tdqcn(i,j,k,2).ge.tdqcn(i,j,k,ntype-1)) ! Si activité de l'espece 2
401     &         tdqcn(i,j,k,ntype-1)=tdqcn(i,j,k,2)
402               if(ival3.eq.1.and.
403     &         tdqcn(i,j,k,3).ge.tdqcn(i,j,k,ntype-1)) ! Si activité de l'espece 3
404     &         tdqcn(i,j,k,ntype-1)=tdqcn(i,j,k,3)
405
406               if(tdqcn(i,j,k,ntype-1).le.-0.99e39)
407     &         tdqcn(i,j,k,ntype-1)=0.
408
409               tdq(i,j,k,ntype-1)=1.e40                    ! plus grande valeur possible
410
411               if(ival1.eq.1 .and. tdq(i,j,k,1).le.tdq(i,j,k,ntype-1)) ! Si activité de l'espece 1
412     &           tdq(i,j,k,ntype-1)=tdq(i,j,k,1)
413               if(ival2.eq.1 .and. tdq(i,j,k,2).le.tdq(i,j,k,ntype-1)) ! Si activité de l'espece 2
414     &           tdq(i,j,k,ntype-1)=tdq(i,j,k,2)
415               if(ival3.eq.1 .and. tdq(i,j,k,3).le.tdq(i,j,k,ntype-1)) ! Si activité de l'espece 3
416     &           tdq(i,j,k,ntype-1)=tdq(i,j,k,3)
417
418               if(tdq(i,j,k,ntype-1).ge.0.99e39) tdq(i,j,k,ntype-1)=0.
419
420                tqcn(i,j,k)=tqcn(i,j,k)+tdqcn(i,j,k,ntype-1)           ! Alors on ajoute les tendances (positive pour qcn ?)
421                tq(i,j,k)  =tq(i,j,k)  +tdq(i,j,k,ntype-1)             !
422
423                if(tqcn(i,j,k).le.0.) tqcn(i,j,k)=0.                   ! et on régularise les tableaux noyaux et aerosols.
424                if(tq(i,j,k)  .le.0.) tq(i,j,k)=0.                     !
425
426             enddo
427           enddo
428         enddo
429
430        continue
431
4321202    format(i2,1x,i2,6(1x,e12.4) )
433       return
434       end
Note: See TracBrowser for help on using the repository browser.