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

Last change on this file since 778 was 175, checked in by slebonnois, 14 years ago

S.LEBONNOIS:

  • Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
  • Quelques ajouts dans la doc.
File size: 13.7 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)         
133         real gg,xmair
134
135         integer jsup,jinf,h,i,j,k,ndim
136         integer ival1,ival2,ival3
137
138         integer iprem
139
140         save iprem,xprime
141         data iprem/0/
142
143         ndim=3*nrad+1
144
145         gg=g0
146
147*********************************************
148*********************************************
149*  Appel de la condensation du methane
150*********************************************
151*********************************************
152
153
154
155*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
156*                   Bilan avant sur le methane                         *
157*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
158
159         do i=1,ng1  ! ng1=1 !!
160           do j=1,NL
161
162*  RAZ des delta d'especes communes.
163*----------------------------------
164             do k=1,nrad
165               tdqcn(i,j,k,1)= 0.
166               tdq( i,j,k,1) = 0.
167             enddo
168
169             gg=g0*rtit**2/(rtit+z(j))**2.
170             xmair=(pb(j+1)-pb(j))/gg/dzb(j)
171
172             do k=1,nrad
173               especes(i,j,k)=tq(i,j,k)         /xmair       ! aerosols, noyaux,
174               especes(i,j,k+nrad)=tqc1(i,j,k)  /xmair       ! methane condense,
175               condens(i,j,k)=(tqc2(i,j,k)+tqc3(i,j,k))/xmair! autre(s) condensat(s)
176               especes(i,j,k+2*nrad)=tqcn(i,j,k)/xmair       ! nombre/kg &  volume/kg
177             enddo
178             especes(i,j,3*nrad+1)=gaz1(i,j)*mch4/mair       ! methane gazeux kg/kg
179           enddo
180         enddo
181
182
1831001   format(7(1x,e12.6),' avN2CH4C2H6')
1841003   format(7(1x,e12.6),' miN2CH4C2H6')
1851002   format(7(1x,e12.6),' apN2CH4C2H6')
186*
187         call n_methane(ng1,ndim,nrad,ddt,
188     &                  p,t,r_e,especes,condens)
189
190
191*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
192*                   Bilan apres sur le methane                         *
193*=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
194
195         do i=1,ng1
196           do j=1,NL
197
198             gg=g0*rtit**2/(rtit+z(j))**2.
199             xmair=(pb(j+1)-pb(j))/gg/dzb(j)
200
201*  ici ce sont les tendances a sortir de nuages.F pour le methane....
202*-------------------------------------------------------------------
203             sum=0.   
204             do k=1,nrad
205               tdqc1(i,j,k)=(especes(i,j,k+nrad)*xmair-tqc1(i,j,k) )
206               sum =sum+tdqc1(i,j,k)/xmair*rhoi_ch4
207               tqc1(i,j,k) = especes(i,j,k+nrad)*xmair
208             enddo
209
210             dgaz1(i,j)= especes(i,j,3*nrad+1)*mair/mch4-gaz1(i,j)
211             gaz1(i,j)= especes(i,j,3*nrad+1)*mair/mch4
212
213c         dgaz1(i,j)=-sum*xmuair/16.
214c          gaz1(i,j)=gaz1(i,j)+dgaz1(i,j)
215
216*  Premiere tendance sur les variables communes (aerosols et noyaux)
217*------------------------------------------------------------------
218             do k=1,nrad
219               tdqcn(i,j,k,1)=(especes(i,j,k+2*nrad)*xmair-tqcn(i,j,k))
220               tdq( i,j,k,1) =(especes(i,j,k)*xmair       -tq(i,j,k)) 
221             enddo
222           enddo
223         enddo
224
225
226
227*     attention, si il y a  de l'ethane sur les noyaux... il est impossible
228*     de les restituer - en revanche on peut en creer de nouveaux ! !!
229*     Le corrolaire de la condition ci dessus est que si il est impossible de
230*     restituer des noyaux, le nombre d'aeorsols ne peux pas augmenter, il
231*     peut en revanche diminuer
232
233
234*********************************************
235*********************************************
236*  Appel de la condensation de l'ethane
237*********************************************
238*********************************************
239
240         do i=1,ng1
241           do j=1,NL
242             do k=1,nrad
243               tdqcn(i,j,k,2)= 0.
244               tdq( i,j,k,2) = 0.
245             enddo
246
247             gg=g0*rtit**2/(rtit+z(j))**2.
248             xmair=(pb(j+1)-pb(j))/gg/dzb(j)
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
269             gg=g0*rtit**2/(rtit+z(j))**2.
270             xmair=(pb(j+1)-pb(j))/gg/dzb(j)
271
272*  ici ce sont les tendances a sortir de nuages.F pour l'ethane....
273*-----------------------------------------------------------------
274
275             sum=0.
276             do k=1,nrad
277               tdqc2(i,j,k)=(especes(i,j,k+nrad)*xmair-tqc2(i,j,k) )
278               sum =sum+tdqc2(i,j,k)/xmair*rhoi_c2h6
279               tqc2(i,j,k) = especes(i,j,k+nrad)*xmair
280             enddo
281
282             dgaz2(i,j)=especes(i,j,3*nrad+1)*mair/mc2h6 - gaz2(i,j)
283             gaz2(i,j)=especes(i,j,3*nrad+1)*mair/mc2h6
284
285c         dgaz2(i,j)=-sum*xmuair/30.
286c          gaz2(i,j)=gaz2(i,j)+dgaz2(i,j)
287
288
289*  Deuxieme tendance sur les variables communes (aerosols et noyaux)
290*------------------------------------------------------------------
291   
292             do k=1,nrad
293               tdqcn(i,j,k,2)=(especes(i,j,k+2*nrad)*xmair-tqcn(i,j,k))
294               tdq(i,j,k,2)  =(especes(i,j,k)*xmair       -tq(i,j,k))
295             enddo
296
297           enddo
298         enddo
299
300
301*********************************************
302*********************************************
303*  Appel de la condensation de l'acethylene
304*********************************************
305*********************************************
306
307         do i=1,ng1
308           do j=1,NL
309             do k=1,nrad
310               tdqcn(i,j,k,3)= 0.
311               tdq( i,j,k,3) = 0.
312             enddo
313             gg=g0*rtit**2/(rtit+z(j))**2.
314             xmair=(pb(j+1)-pb(j))/gg/dzb(j)
315
316             do k=1,nrad
317               especes(i,j,k)=tq(i,j,k)/xmair
318               especes(i,j,k+nrad)=tqc3(i,j,k)/xmair           ! acethylene condense
319               condens(i,j,k)=(tqc1(i,j,k)+tqc2(i,j,k))/xmair  ! autres condensats
320               especes(i,j,k+2*nrad)=tqcn(i,j,k)/xmair
321             enddo
322
323             especes(i,j,3*nrad+1)=gaz3(i,j)*mc2h2/mair      ! acethylene gazeux
324
325           enddo
326         enddo
327
328         call n_acethylene(ng1,ndim,nrad,ddt,
329     &                     p,t,r_e,especes,condens)
330
331
332         do i=1,ng1
333           do j=1,NL
334
335             gg=g0*rtit**2/(rtit+z(j))**2.
336             xmair=(pb(j+1)-pb(j))/gg/dzb(j)
337
338*  ici ce sont les tendances a sortir de nuages.F pour l'ethane....
339*-----------------------------------------------------------------
340
341             sum=0.
342             do k=1,nrad
343               tdqc3(i,j,k)=(especes(i,j,k+nrad)*xmair-tqc3(i,j,k) )
344               sum =sum+tdqc3(i,j,k)/xmair*rhoi_c2h2
345               tqc3(i,j,k) = especes(i,j,k+nrad)*xmair
346             enddo
347
348             dgaz3(i,j)=especes(i,j,3*nrad+1)*mair/mc2h2 - gaz3(i,j)
349             gaz3(i,j)=especes(i,j,3*nrad+1)*mair/mc2h2
350
351c         dgaz3(i,j)=-sum*xmuair/26.
352c         gaz3(i,j)=gaz3(i,j)+dgaz3(i,j)
353
354*  Troisieme tendance sur les variables communes (aerosols et noyaux)
355*------------------------------------------------------------------
356
357
358             do k=1,nrad
359               tdqcn(i,j,k,3)=(especes(i,j,k+2*nrad)*xmair-tqcn(i,j,k))
360               tdq(i,j,k,3)  =(especes(i,j,k)*xmair       -tq(i,j,k))
361             enddo
362
363           enddo
364         enddo
365
366
367*  FIN DES APPELS DE NUAGES ET BILAN DES TENDANCES...
368*------------------------------------------------------------------
369
370
371         do i=1,ng1
372           do j=1,NL
373             do k=1,nrad
374
375*          Ici on test l'activité nuageuse : si on a l'association
376*          tdqcX(i,j,k) = 0 et tqcX(i,j,k) = 0 alors ivalX reste à 0 (pas d'ativité)
377*                                              sinon ivalX passe à 1 (activité)
378*------------------------------------------------------------------------------------
379
380               ival1=0
381               ival2=0
382               ival3=0
383
384               if(tdqc1(i,j,k).ne.0. .or. tqc1(i,j,k).gt.0.)  ival1=1
385               if(tdqc2(i,j,k).ne.0. .or. tqc2(i,j,k).gt.0.)  ival2=1
386               if(tdqc3(i,j,k).ne.0. .or. tqc3(i,j,k).gt.0.)  ival3=1
387
388*          Ici on definit la tendances des noyaux en faisant deux choses:
389*          -1 On ecarte les cas  tdqcn(i,j,k,X)=0 si ils sont associés à une
390*          absence d'activité nuageuse de l'espèce (tdqcX(i,j,k)=0.)
391*          -2 Sélectionne la tendance la plus élevée. Si aucune activité nuageuse
392*           n'exits dans cette case (ivalX=0 pour les 3 especes), alors on
393*           retrouve la valeur -1.e40 que l'on mets alors à 0.
394*----------------------------------------------------------------------
395c23456789012345678901234567890123456789012345678901234567890123456789012
396
397               tdqcn(i,j,k,ntype-1)=-1.e40                      ! plus petite valeur possible
398
399               if(ival1.eq.1.and.
400     &         tdqcn(i,j,k,1).ge.tdqcn(i,j,k,ntype-1)) ! Si activité de l'espece 1
401     &         tdqcn(i,j,k,ntype-1)=tdqcn(i,j,k,1)
402               if(ival2.eq.1.and.
403     &         tdqcn(i,j,k,2).ge.tdqcn(i,j,k,ntype-1)) ! Si activité de l'espece 2
404     &         tdqcn(i,j,k,ntype-1)=tdqcn(i,j,k,2)
405               if(ival3.eq.1.and.
406     &         tdqcn(i,j,k,3).ge.tdqcn(i,j,k,ntype-1)) ! Si activité de l'espece 3
407     &         tdqcn(i,j,k,ntype-1)=tdqcn(i,j,k,3)
408
409               if(tdqcn(i,j,k,ntype-1).le.-0.99e39)
410     &         tdqcn(i,j,k,ntype-1)=0.
411
412               tdq(i,j,k,ntype-1)=1.e40                    ! plus grande valeur possible
413
414               if(ival1.eq.1 .and. tdq(i,j,k,1).le.tdq(i,j,k,ntype-1)) ! Si activité de l'espece 1
415     &           tdq(i,j,k,ntype-1)=tdq(i,j,k,1)
416               if(ival2.eq.1 .and. tdq(i,j,k,2).le.tdq(i,j,k,ntype-1)) ! Si activité de l'espece 2
417     &           tdq(i,j,k,ntype-1)=tdq(i,j,k,2)
418               if(ival3.eq.1 .and. tdq(i,j,k,3).le.tdq(i,j,k,ntype-1)) ! Si activité de l'espece 3
419     &           tdq(i,j,k,ntype-1)=tdq(i,j,k,3)
420
421               if(tdq(i,j,k,ntype-1).ge.0.99e39) tdq(i,j,k,ntype-1)=0.
422
423                tqcn(i,j,k)=tqcn(i,j,k)+tdqcn(i,j,k,ntype-1)           ! Alors on ajoute les tendances (positive pour qcn ?)
424                tq(i,j,k)  =tq(i,j,k)  +tdq(i,j,k,ntype-1)             !
425
426                if(tqcn(i,j,k).le.0.) tqcn(i,j,k)=0.                   ! et on régularise les tableaux noyaux et aerosols.
427                if(tq(i,j,k)  .le.0.) tq(i,j,k)=0.                     !
428
429             enddo
430           enddo
431         enddo
432
433        continue
434
4351202    format(i2,1x,i2,6(1x,e12.4) )
436       return
437       end
Note: See TracBrowser for help on using the repository browser.