source: trunk/LMDZ.TITAN/libf/phytitan/n_methane.F @ 201

Last change on this file since 201 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: 16.2 KB
Line 
1      subroutine n_methane(ngrid,nq,nbin,
2     *                     dt,pl,tl,aerad,
3     *                     q,qprime)
4
5      implicit none
6#include "dimensions.h"
7#include "microtab.h"
8#include "varmuphy.h"
9
10c  Arguments
11c  ---------
12
13      integer ngrid,nq,nbin
14
15      REAL dt              ! physical time step (s)
16      REAL pl(ngrid,nz)    ! pressure at each level (mbar)
17      REAL tl(ngrid,nz)    ! temperature at each level (K)
18      REAL aerad(nbin)     ! Radius array
19
20c    Tracers :
21      REAL q(ngrid,nz,nq)         ! tracer (kg/kg)
22      REAL qprime(ngrid,nz,nbin)  ! tracer (kg/kg)
23
24
25c     Local variables
26c     ---------------
27
28      integer ntyp
29      parameter (ntyp=3)   ! 3 = 1 aerosol + 1 noyau + 1 glace
30
31      real n_aer(nz,nbin,ntyp)
32      real ch4vap(nz)
33
34      integer itrac
35      integer ig,i,j,k,l,n   ! Loop integers
36      integer ilay,iq
37
38c  Treatment
39c  ---------
40
41      DO ig = 1 , NGRID
42
43c     Set up the aerosol array
44        do j = 1, ntyp
45          do k = 1, nbin
46            itrac = (j-1) * nbin + k
47            do l = 1, nz
48              n_aer(l,k,j) = max(q(ig,l,itrac),0.)
49            enddo
50          enddo
51        enddo
52
53c     Set up the methane vapor array
54        do l = 1, nz
55          ch4vap(l) = q(ig,l,nq)
56        enddo
57
58
59
60        call nucleacond(ngrid,nbin,dt,ig,pl,tl,aerad,
61     &                 n_aer,qprime,ch4vap)
62
63
64c       Update q arrays
65        do j = 1, ntyp
66          do k = 1, nbin
67            itrac = (j-1) * nbin + k
68            do l = 1, nz
69              q(ig,l,itrac) = n_aer(l,k,j)
70            enddo
71          enddo
72        enddo
73
74c     Update methane vapor array
75        do l = 1, nz
76          q(ig,l,nq) = ch4vap(l)
77        enddo
78
79      ENDDO
80
81      return
82      END
83
84****************************************************************
85      subroutine nucleacond(ngrid,nbin,dt,ig,
86     &                      pl,tl,aerad,n_aer,qprime,ch4vap)
87*                                                              *
88*     This routine updates species concentrations due          *
89*     to both nucleation and condensation-induced variations.  *
90*     Gain and loss rates associated to each one of these      *
91*     processes are computed separately in other routines.     *
92*                                                              *
93****************************************************************
94
95      implicit none
96#include "dimensions.h"
97#include "microtab.h"
98#include "varmuphy.h"
99
100      integer ng,nalt
101      parameter(ng=1,nalt=llm)
102
103      real lv
104      common/lheat/lv
105
106c  Arguments
107c  ---------
108
109      integer ngrid,nbin
110      integer ig
111      integer ntyp
112      parameter (ntyp=3)
113
114      real dt                    ! Global time step
115      real pl(ngrid,nz),tl(ngrid,nz)
116      real aerad(nbin)
117      real ch4vap(nz)            ! Methane vapor mass mixing ratio (kg/m3)
118      real ch4vap_old
119      real n_aer(nz,nbin,ntyp)  ! number concentrations of particle/each size bin
120      real qprime(ngrid,nz,nbin)  ! tracer (kg/kg)
121      REAL total1(nz),total11(nz),total2(nz),total22(nz)
122      REAL dmsm,mtot
123
124
125c  Local
126c  -----
127
128      integer i,j,k,l,n,iindice,iselec
129
130      real dQc           ! Amount of condensed methane (kg/m3) during timestep
131      real*8 sat_ratio   ! Methane saturation ratio over liquid
132      real*8 sat_ratmix   ! Methane saturation ratio over liquid
133      real*8 pch4        ! Methane partial pressure (Pa)
134      real qsat          ! Methane mass mixing ratio at saturation (kg/kg of air)
135      real qsatmix          ! Methane mass mixing ratio at saturation (kg/kg of air)
136      real*8 rate(nbin)    ! Heterogeneous Nucleation rate (s-1)
137      real*8 elim         
138
139      real nsav(nbin,ntyp)
140      real dn(nbin,ntyp)
141      real rad(nbin)     ! Radius of droplets in each size bin
142      real*8 gr(nbin)      ! Growth rate in each bin
143      real radius        ! Radius of droplets after growth
144      real Qs            ! Mass of condensate required to reach saturation
145      real newsat
146      real vol(nbin)
147
148      real press
149      real sig,temp,seq(nbin)
150      real Ctot,up,dwn,newvap,gltot
151
152      real rho_a,cap
153      real tempref
154      real frac
155      real xtime,xtime_prime
156      real dQc2
157
158c     Variables for latent heat release
159      real lw,cpp
160      data lw / 510.e+3/
161c     data cpp/1044./
162      data cpp/1050./ ! pour etre cohérent avec le reste...
163      save lw,cpp
164
165
166c  Treatment
167c  ---------
168
169c  Treatment
170c  ---------
171      do i = 1, nbin
172        vol(i) = 4./3. * pi * aerad(i)**3.
173      enddo
174
175        do l = 1, nz
176        total1(l)=0. !solide
177        do k = 1, nbin
178        total1(l)=total1(l)+n_aer(l,k,2)*rhoi_ch4
179        enddo
180        total2(l)=ch4vap(l)
181        enddo
182
183c     Start loop over heights
184      DO 100 l = 1, nz
185   
186        iindice=0                ! mettre l'indice à 0
187
188        temp   = tl(ig,l)
189        press  = pl(ig,l)
190        tempref=temp
191
192c       Save the values of the particle arrays before condensation
193        do j = 1, ntyp
194          do i = 1, nbin
195            nsav(i,j) = n_aer(l,i,j)
196          enddo
197        enddo
198
199
200 99     continue   ! DEBUT DE LA BOUCLE CONDITONNELLE
201
202        call ch4sat(temp,press,qsat)
203        qsat=qsat*0.85
204        qsatmix=qsat*0.85
205
206c       Get the partial presure of methane vapor and its saturation ratio
207        pch4       = ch4vap(l) * (Mn2/Mch4) * press
208        sat_ratio  = ch4vap(l) / qsat
209        sat_ratmix = ch4vap(l) / qsatmix
210
211c       Get the rates of nucleation
212        call nuclea(nbin,aerad,pch4,temp,sat_ratio,rate)
213
214c       Get the growth rates of condensation/sublimation
215        up   = ch4vap(l)
216        dwn  = 1.
217        Ctot = ch4vap(l)
218        DO i = 1, nbin
219
220        if (n_aer(l,i,3).eq.0) then
221          rad(i) = aerad(i)
222        else
223          rad(i) = ((n_aer(l,i,2)/n_aer(l,i,3) +
224     &     qprime(ig,l,i)/n_aer(l,i,3)
225     &    +vol(i))*0.75/pi)**(1./3.)
226        endif
227       
228
229*       Equilibrium saturation ratio (due to curvature effect)
230        seq(i) = exp( 2.*sig(temp)*Mch4 / (rhoi_ch4*rgp*temp*rad(i)) )
231
232        call growthrate(dt,temp,press,pch4,
233     &                  sat_ratmix,seq(i),rad(i),gr(i))
234        up = up + dt * gr(i) * 4. * pi * rhoi_ch4 * rad(i) * seq(i)
235     *                 * nsav(i,3)
236        dwn= dwn+ dt * gr(i) * 4. * pi * rhoi_ch4 * rad(i) / qsat
237     *                 * nsav(i,3)
238        Ctot= Ctot + rhoi_ch4 * nsav(i,2)
239
240        ENDDO
241
242        newvap = min(up/dwn,Ctot) 
243        newvap = max(newvap,0.) 
244
245        gltot=0.
246        DO i = 1, nbin
247          gr(i)  = gr(i) * ( newvap/qsat - seq(i) )
248          if(nsav(i,2).le.0. .and. gr(i).le.0.) then
249              n_aer(l,i,2) = 0.
250          else
251          n_aer(l,i,2) = nsav(i,2) + dt * gr(i) * 4. * pi * rad(i)
252     *                               * n_aer(l,i,3)
253          if (n_aer(l,i,2).le.0.) then
254            n_aer(l,i,1) = n_aer(l,i,1) + n_aer(l,i,3)
255            n_aer(l,i,2) = 0.
256            n_aer(l,i,3) = 0.
257          endif
258          gltot=n_aer(l,i,2)*rhoi_ch4+gltot
259          endif
260
261        ENDDO
262   
263c       Determine the mass of exchanged methane
264
265        dQc = 0.
266        DO i = 1, nbin
267          dQc = dQc - rhoi_ch4 * ( n_aer(l,i,2) - nsav(i,2) )
268        ENDDO
269
270
271c       Arrays resetted to their initial value before condensation
272
273        do j = 1, ntyp
274          do i = 1, nbin
275            dn(i,j)      = n_aer(l,i,j) - nsav(i,j)
276            n_aer(l,i,j) = nsav(i,j)
277          enddo
278        enddo
279
280
281c       Update the c arrays.
282c       nucleation & cond/eva tendencies added together.
283
284        do i=1,nbin
285          elim         = dt * rate(i)
286          n_aer(l,i,1) = n_aer(l,i,1) / (1.+elim)
287          n_aer(l,i,3) = n_aer(l,i,3) + elim * n_aer(l,i,1) + dn(i,3)
288          n_aer(l,i,1) = n_aer(l,i,1) + dn(i,1)
289          n_aer(l,i,2) = n_aer(l,i,2) + dn(i,2)
290         if(n_aer(l,i,2).lt.0.) n_aer(l,i,2)=0.
291
292        enddo
293
294        dQc = 0.
295        DO i = 1, nbin    ! dQc <0 si glace produite !!
296          dQc = dQc - rhoi_ch4 * ( n_aer(l,i,2) - nsav(i,2) )
297        ENDDO
298
299
300       ch4vap(l)    = ch4vap(l) + dQc
301
302100   CONTINUE
303
304        do l = 1, nz
305        total11(l)=0.
306        do k = 1, nbin
307          total11(l)=total11(l)+n_aer(l,k,2)*rhoi_ch4
308        enddo
309        total22(l)=ch4vap(l)
310        enddo
311
312
313      return
314      end
315
316
317*******************************************************
318*                                                     *
319      subroutine nuclea(nbin,aerad,pch4,temp,sat,nucrate)
320*   This subroutine computes the nucleation rate      *
321*   as given in Pruppacher & Klett (1978) in the      *
322*   case of water ice forming on a solid substrate.   *
323*     Definition refined by Keese (jgr,1989)          *
324*                                                     *
325*******************************************************
326
327      implicit none
328#include "dimensions.h"
329#include "microtab.h"
330#include "varmuphy.h"
331
332      integer nbin
333      real aerad(nbin)
334
335      real*8 nucrate(nbin)
336      real*8 pch4
337      real   temp
338      real*8 sat
339
340      integer l,i
341      real*8 nch4
342      real sig            ! Water-ice/air surface tension  (N.m)
343      real*8 rstar        ! Radius of the critical germ (m)
344      real*8 gstar        ! # of molecules forming a critical embryo
345      real*8 x            ! Ratio rstar/radius of the nucleating dust particle
346      real fistar         ! Activation energy required to form a critical embryo (J)
347      real*8 zeldov       ! Zeldovitch factor (no dim)
348      real*8 fshape       ! function defined at the end of the file
349      real*8 deltaf
350
351      real nus
352      data nus/1.e+13/       ! Jump frequency of a molecule (s-1)
353      real m0
354      data m0/2.6578e-26/     ! Weight of a methane molecule (kg)
355      real vo1
356c     data vo1/6.254e-29/     ! Volume of a methane molecule (m3)
357      data vo1/3.7649e-5/     ! Volume of a methane molecule (m3)
358
359      real desorp
360      data desorp/0.288e-19/ ! Activation energy for desorption of water on a dust-like substrate (J/molecule)
361      real surfdif
362      data surfdif/0.288e-20/! Estimated activation energy for surface diffusion of water molecules (J/molecule)
363
364      IF (sat .GT. 1.) THEN    ! minimum condition to activate nucleation
365
366        nch4    = pch4 / kbz / temp
367        rstar  = 2. * sig(temp) * vo1 / (rgp*temp*log(sat))
368        gstar  = 4. * nav * pi * (rstar**3) / (3.*vo1)
369c       Loop over size bins
370        do i=1,nbin
371          x      = aerad(i) / rstar
372          x      = aerad(imono) / rstar  ! attention r(5)=monomeres
373
374          fistar = (4./3.*pi) * sig(temp) * (rstar**2.) *
375     &             fshape(mtetach4,x)
376          deltaf = min( max((2.*desorp-surfdif-fistar)/(kbz*temp)
377     &           , -100.), 100.)
378          if (deltaf.eq.-100.) then
379            nucrate(i) = 0.
380          else
381            zeldov = sqrt ( fistar / (3.*pi*kbz*temp*(gstar**2.)) )
382            nucrate(i)  = zeldov * kbz* temp * rstar**2.
383     &                  * 4. * pi * ( nch4*aerad(i) )**2.
384     &                  / ( fshape(mtetach4,x) * nus * m0 )
385     &                  * dexp(deltaf)
386 
387             if(i.gt.imono)  nucrate(i)= zeldov * kbz* temp * rstar**2.
388     &           * 4. * pi * vrat_e**(i-imono)*(nch4*aerad(imono) )**2.
389     &            / (fshape(mtetach4,x) * nus * m0 )
390     &            * dexp(deltaf)
391
392
393          endif
394        enddo
395      ELSE
396        do i=1,nbin
397          nucrate(i) = 0.
398        enddo
399
400      ENDIF
401
402      return
403      end
404
405******************************************************************
406        subroutine ch4sat(t,p,qsat)
407* cette fonction calcule la pression de vapeur saturante du      *
408* methane a une altitude donnee z par l'equation de Thek-Stiel   *
409******************************************************************
410
411        real rgp
412        data rgp/8.3143/
413
414* declaration des variables internes
415* ----------------------------------
416        real p,tc,pc,tr,tb,tbr,phib,ac,hbr,hvb,avb,a,b
417        real qsat
418
419        tc = 190.53
420        pc = 45.96 * 0.986923
421        tr = t / tc
422        tb = 111.63
423        tbr= tb / tc
424        phib = -35. + 36./tbr + 42.*log(tbr)-tbr**6
425        ac   = (0.315*phib + log(pc))/(0.0838*phib-log(tbr))
426        hbr  = tbr * log(pc) / (1.-tbr)
427        hvb  = rgp*tb*(0.4343*log(pc)-0.68859+0.89584*tbr)/(0.37691-
428     s      0.37306*tbr+0.14878/(pc*tbr**2))
429        avb  = hvb / (rgp*tc*(1.-tbr)**(0.375))
430        a    = 5.2691 + 2.0753*avb - 3.1738*hbr
431        b    = 1.042 * ac - 0.46284 * avb
432
433        qsat=pc*exp(avb*(1.14893-1./tr-0.11719*tr-0.03174*tr**2
434     s        -0.375*log(tr))+b*((tr**a-1.)/a+0.04*(1./tr-1.)))
435        qsat=qsat*1e+5/0.986923
436        !  ---> qsat en Pa
437
438        qsat=qsat* 16.0 / (28.0*p)   ! kg/kg
439       
440
441
442        return
443        end
444
445c=======================================================================
446      subroutine growthrate(timestep,temp,press,pch4,sat,seq,r,Cste)
447c
448c     Determination of the droplet growth rate
449c
450c=======================================================================
451
452      IMPLICIT NONE
453#include "dimensions.h"
454#include "microtab.h"
455#include "varmuphy.h"
456
457c-----------------------------------------------------------------------
458c   declarations:
459c   -------------
460
461      common/lheat/Lv
462
463c
464c   arguments:
465c   ----------
466
467      REAL timestep
468      REAL temp    ! temperature in the middle of the layer (K)
469      REAL press    ! pressure in the middle of the layer (K)
470      REAL*8 pch4 ! Methane vapor partial pressure (Pa)
471      REAL*8 sat  ! saturation ratio
472      REAL r    ! crystal radius before condensation (m)
473      REAL seq  ! Equilibrium saturation ratio
474
475c   local:
476c   ------
477
478      REAL psat
479      REAL moln2,molch4
480      REAL To,wch4,tch4,ftm
481
482c     Effective gas molecular radius (m)
483      data moln2/1.75e-10/   ! N2
484c     Effective gas molecular radius (m)
485      data molch4/2.e-10/   ! CH4
486
487      data tch4/190.4/
488      data wch4/1.1e-2/    ! Reid et al (Eq 7-9.4 + Appendix Compound [116])
489
490      REAL k,Lv                 
491      REAL knudsen           ! Knudsen number (gas mean free path/particle radius)
492      REAL a,Dv,lambda,Rk,Rd ! Intermediate computations for growth rate
493      REAL*8 Cste
494
495c-----------------------------------------------------------------------
496c      Ice particle growth rate by diffusion/impegement of molecules
497c                r.dr/dt = (S-Seq) / (Seq*Rk+Rd)
498c        with r the crystal radius, Rk and Rd the resistances due to
499c        latent heat release and to vapor diffusion respectively
500c-----------------------------------------------------------------------
501
502      psat = pch4 / sat
503
504c     - Thermal conductibility of N2
505      k = ( 2.857e-2 * temp - 0.5428  ) * 4.184e-3
506c     - Latent heat of ch4 (J.kg-1)
507      Lv = 510.e+3
508      ftm=1.-temp/tch4
509      if (ftm.le.1.e-3)  ftm=1.e-3
510      Lv=8.314*tch4*(7.08*ftm**0.354+10.95*wch4*ftm**0.456)/16.e-3
511     
512
513c     - Constant to compute gas mean free path
514c     l= (T/P)*a, with a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
515
516      a = 0.707*rgp/(4 * pi* (moln2*1.e10)**2  * (nav*1.e-20))
517
518c     - Compute Dv, methane vapor diffusion coefficient
519c       accounting for both kinetic and continuum regime of diffusion,
520c       the nature of which depending on the Knudsen number.
521
522      Dv = 1./3. * sqrt( 8*rgp*temp/(pi*Mch4) )* (kbz*1.e20) * temp /
523     &   (pi*press*(moln2*1.e10+molch4*1.e10)**2 * sqrt(1.+Mch4/Mn2) )
524
525      knudsen = temp / press * a / r
526      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen)
527      Dv      = Dv / (1. + lambda * knudsen)
528
529c     - Compute Rk
530      Rk = Lv**2 * rhoi_ch4 * Mch4 / (k*rgp*temp**2.)
531c     - Compute Rd
532      Rd = rgp * temp *rhoi_ch4 / (Dv*psat*Mch4)
533
534c     - Compute:      rdr/dt = Cste * (S-Seq)
535      Cste = 1. / (seq*Rk+Rd)
536
537      RETURN
538      END
539
540
541*********************************************************
542      real function sig(t)
543*    this function computes the surface tension (N.m)   *
544*   between methane and air as a function of temp.    *
545*********************************************************
546
547      real t
548
549      sig = ( t/4. + 41.) * 1e-3
550
551      return
552      end
553
554*********************************************************
555      real*8 function fshape(cost,rap)
556*        function computing the f(m,x) factor           *
557* related to energy required to form a critical embryo  *
558*********************************************************
559
560      implicit none
561
562      real cost
563      real*8 rap
564      real*8 phi
565      real*8 a,b,c
566
567      phi = sqrt( 1. - 2.*cost*rap + rap**2 )
568
569      a = 1. + ( (1.-cost*rap)/phi )**3
570      b = (rap**3) * (2.-3.*(rap-cost)/phi+((rap-cost)/phi)**3)
571      c = 3. * cost * (rap**2) * ((rap-cost)/phi-1.)
572
573      fshape = 0.5*(a+b+c)
574
575      if (rap.gt.3000.) fshape = ((2.+cost)*(1.-cost)**2)/4.
576
577
578      return
579      end
580
Note: See TracBrowser for help on using the repository browser.