source: trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F @ 1652

Last change on this file since 1652 was 1620, checked in by emillour, 9 years ago

Mars GCM:
Small fix for gfortran, more strict than ifort: logicals should be compared using .eqv., not .eq.
EM
  • Property svn:executable set to *
File size: 17.2 KB
Line 
1
2
3c=======================================================================
4      subroutine massflowrateCO2(P,T,Sat,Radius,Matm,Ic)
5c
6c     Determination of the mass transfer rate
7c
8c     newton-raphson method
9
10c     CLASSICAL  (no SF etc.)
11     
12c     AUTOMATIC SETTING OF RANGES FOR NEWTON-RAPHSON FOR THE PAPER
13
14c     MASS FLUX Ic
15
16c=======================================================================
17      USE comcstfi_h
18
19      implicit none
20
21      include "microphys.h"
22c      include "microphysCO2.h"
23
24
25
26c   arguments: INPUT
27c   ----------
28
29      REAL T,Matm
30      REAL*8 SAT
31      real P
32      DOUBLE PRECISION Radius
33c   arguments: OUTPUT
34c   ----------
35
36      DOUBLE PRECISION   Ic
37
38c   Local Variables
39c   ----------
40
41      DOUBLE PRECISION   Tcm
42      DOUBLE PRECISION   T_inf, T_sup, T_dT
43      DOUBLE PRECISION   C0,C1,C2
44      DOUBLE PRECISION   kmix,Lsub,cond
45      DOUBLE PRECISION   rtsafe
46      DOUBLE PRECISION   left, fval, dfval
47
48c    function for newton-raphson iterative method
49c    --------------------------
50
51      EXTERNAL classical
52
53         
54      Tcm = dble(T)    ! initialize pourquoi 0 et pas t(i)
55
56      T_inf = 0d0
57      T_sup = 200d0
58
59      T_dT = 0.1  ! precision - mettre petit et limiter nb iteration?
60     
61666   CONTINUE
62
63c      print*, 'Radius ', Radius
64c      print*, 'SAT = ', Sat
65      call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2)
66
67      if (isnan(C0) .eqv. .true.)  C0=0d0
68
69   
70c     FIND SURFACE TEMPERATURE (Tc) : iteration sur t
71
72      cond = 4.*pi*Radius*kmix
73 
74      Tcm = rtsafe(classical,T_inf,T_sup,T_dT,Radius,C0,C1,C2)
75
76
77      if (Tcm.LE.0d0) then ! unsignificant cases where S<<<Seq and Ncores <<1e-10
78
79            Tcm = 0d0
80
81      endif
82
83
84c     THEN COMPUTE MASS FLUX Ic from FINAL Tsurface (Tc)
85
86      Ic = (Tcm-T)
87
88      Ic = cond*Ic/(-Lsub)
89c regarder de combien varie la solution Ic entre Tcm et Tcm+T_dT
90   
91      RETURN
92
93      END
94
95
96c****************************************************************
97
98      FUNCTION rtsafe(funcd,x1,x2,xacc,Radius,C0,C1,C2)
99*
100*
101*     Newton Raphsen routine (Numerical Recipe)
102*
103c****************************************************************
104
105      implicit none
106
107      INTEGER MAXIT
108       DOUBLE PRECISION x1,x2,xacc
109      DOUBLE PRECISION rtsafe
110      DOUBLE PRECISION Radius
111      DOUBLE PRECISION C0,C1,C2
112
113
114      EXTERNAL funcd
115
116      PARAMETER (MAXIT=10000)   !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection,
117                                !find the root of a function bracketed between x1 and x2. The root, returned as the function value rtsafe,
118                                !will be refined until its accuracy is known within !!±xacc. funcd is a user-supplied subroutine which
119                                !returns both the function value and the first derivative of the function.
120
121      INTEGER j
122      DOUBLE PRECISION df,dx,dxold
123      DOUBLE PRECISION f,fh,fl,temp,xh,xl
124
125      call funcd(x1,fl,df,C0,C1,C2)
126      call funcd(x2,fh,df,C0,C1,C2)
127
128
129      if ((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.) ) then
130         
131         
132         x1=0d0
133         x2=500d0
134         
135         call funcd(x1,fl,df,C0,C1,C2)
136         call funcd(x2,fh,df,C0,C1,C2)
137         
138         write(*,*) 'root must be bracketed in rtsafe'
139      endif
140
141
142      if (fl.eq.0.) then
143
144         rtsafe=x1
145         return
146
147      else if (fh.eq.0.) then
148
149         rtsafe=x2
150         return
151
152      else if (fl.lt.0.) then   !Orient the search so that f(xl) < 0.
153                                                                                               
154         xl=x1
155         xh=x2
156
157      else
158
159        xh=x1
160        xl=x2
161
162      endif
163
164      rtsafe = .5*(x1+x2)       !Initialize the guess for root,
165      dxold  = abs(x2-x1)       !the stepsize before last,
166      dx     = dxold            ! and the last step.
167
168
169      call funcd(rtsafe,f,df,C0,C1,C2)
170
171      DO 11 j=1,MAXIT           !Loop over allowed iterations.
172
173         
174         !print*, 'iteration:', j
175         !print*, rtsafe
176
177
178         if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).gt.0.   ! Bisect if Newton out of range
179     *    .or. abs(2.*f).gt.abs(dxold*df) ) then               ! or not decreasing fst enough
180
181             dxold=dx
182             dx=0.5*(xh-xl)
183             rtsafe=xl+dx
184
185             if (xl.eq.rtsafe) return                   !Change in root is negligible. Newton step acceptable. Take it.
186
187         else
188
189            dxold=dx
190            dx=f/df
191            temp=rtsafe
192
193            rtsafe=rtsafe-dx
194
195            if(temp.eq.rtsafe) return
196
197         endif
198
199
200        if(abs(dx).lt.xacc) return !Convergence criterion. The one new function evaluation per iteration. Maintain the bracket on the root.
201
202         call funcd(rtsafe,f,df,C0,C1,C2)
203
204        if(f.lt.0.) then
205         xl=rtsafe
206        else
207         xh=rtsafe
208        endif
209
210
21111    ENDDO
212
213      write(*,*) 'rtsafe exceeding maximum iterations'
214
215      return
216
217      END
218
219
220c********************************************************************************
221
222      subroutine classical(x,f,df,C0,C1,C2)
223             
224c     Function to give as input to RTSAFE (NEWTON-RAPHOEN)
225
226
227c********************************************************************************
228
229      implicit none
230
231      DOUBLE PRECISION   x
232      DOUBLE PRECISION  C0,C1,C2
233
234      DOUBLE PRECISION f
235      DOUBLE PRECISION  df
236     
237
238     
239
240      f   = x + C0*exp(C1*x) - C2   ! start f
241      df  = 1. + C0*C1*exp(C1*x)    ! start df
242 
243      return
244
245      END 
246
247c********************************************************************************
248
249      subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2)
250
251c********************************************************************************
252c defini la fonction eq 6 papier 2014
253      use tracer_mod, only: rho_ice_co2
254      USE comcstfi_h
255
256      implicit none
257
258      include "microphys.h"
259c      include "microphysCO2.h"
260   
261
262c   arguments: INPUT
263c   ----------------
264      REAL P
265      real T
266      REAL*8 S
267      double precision rc
268      REAL Matm !g.mol-1 ( = mmean(ig,l) )
269
270c   local:
271c   ------
272
273
274      DOUBLE PRECISION Cpatm,Cpn2,Cpco2
275      DOUBLE PRECISION psat, xinf, pco2
276      DOUBLE PRECISION Dv     
277      DOUBLE PRECISION l0,l1,l2,l3,l4           
278      DOUBLE PRECISION knudsen, a, lambda      ! F and S correction
279      DOUBLE PRECISION Ak                       ! kelvin factor   
280      DOUBLE PRECISION vthatm,lpmt,rhoatm, vthco2 ! for Kn,th
281
282c   arguments: OUTPUT
283c   ----------
284      DOUBLE PRECISION  C0,C1,C2
285      DOUBLE PRECISION  kmix,Lsub
286
287c     DEFINE heat cap. J.kg-1.K-1 and To
288
289      data Cpco2/0.7e3/
290      data Cpn2/1e3/
291
292      kmix = 0d0
293      Lsub = 0d0
294
295      C0 = 0d0
296      C1 = 0d0
297      C2 = 0d0
298
299c     Equilibirum pressure over a flat surface
300
301      psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T))  ! (Pa)
302
303c     Compute transport coefficient
304
305      pco2 = psat * dble(S)
306
307c     Latent heat of sublimation if CO2  co2 (J.kg-1)
308c     version Azreg_Ainou (J/kg) :
309
310      l0=595594.     
311      l1=903.111     
312      l2=-11.5959   
313      l3=0.0528288
314      l4=-0.000103183
315 
316      Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 *
317     &     dble(T)**3 + l4 * dble(T)**4 ! J/kg
318
319c     atmospheric density
320
321      rhoatm = dble(P*Matm)/(rgp*dble(T))   ! g.m-3
322      rhoatm = rhoatm * 1.00e-3 !kg.m-3
323
324      call  KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2
325      call  Diffcoeff(P, T, Dv) 
326
327      Dv = Dv * 1.00e-4         !!! cm2.s-1  to m2.s-1
328
329c     ----- FS correction for Diff
330
331      vthco2  = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1
332
333      knudsen = 3*Dv / (vthco2 * rc)
334
335      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
336     
337      Dv      = Dv / (1. + lambda * knudsen)
338       
339c     ----- FS correction for Kth
340
341      vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg
342
343      Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1
344
345      lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/
346     &     (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer
347
348      knudsen = lpmt / rc
349
350      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
351
352      kmix    = kmix /  (1. + lambda * knudsen)
353
354c     --------------------- ASSIGN coeff values for FUNCTION
355
356      xinf = dble(S) * psat / dble(P)
357
358      Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) ))
359 
360      C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/
361     &     (rgp*dble(T)))
362 
363      C1 = Lsub*mco2/(rgp*dble(T)**2)
364
365      C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T))
366     
367      RETURN
368      END
369
370
371c======================================================================
372
373      subroutine Diffcoeff(P, T, Diff)
374
375c     Compute diffusion coefficient CO2/N2
376c     cited in Ilona's lecture - from Reid et al. 1987
377c======================================================================
378
379
380       IMPLICIT NONE
381
382       include "microphys.h"
383c       include "microphysCO2.h"
384
385c      arguments
386c     -----------
387     
388      REAL P
389      REAL Pbar                     !!! has to be in bar for the formula
390      REAL T
391     
392c     output
393c     -----------
394
395      DOUBLE PRECISION Diff
396   
397c      local
398c     -----------
399
400      DOUBLE PRECISION dva, dvb, Mab  ! Mab has to be in g.mol-1
401       
402        Pbar = P * 1d-5
403   
404        Mab = 2. / ( 1./mn2 + 1./mco2 ) * 1000.
405
406        dva = 26.9        ! diffusion volume of CO2,  Reid et al. 1987 (cited in Ilona's lecture)
407        dvb = 18.5        ! diffusion volume of N2
408   
409        Diff  = 0.00143 * dble(T)**(1.75) / (dble(Pbar) * sqrt(Mab)
410     &       * (dble(dva)**(1./3.) + dble(dvb)**(1./3.))**2.) !!! in cm2.s-1
411 
412       RETURN
413
414       END
415
416
417c======================================================================
418
419         subroutine KthMixNEW(Kthmix,T,x,rho)
420
421c        Compute thermal conductivity of CO2/N2 mixture
422c         (***WITHOUT*** USE OF VISCOSITY)
423
424c          (Mason & Saxena, 1958 - Wassiljeva 1904)
425c======================================================================
426
427
428       implicit none
429       
430       include "microphys.h"
431c       include "microphysCO2.h"
432
433c      arguments
434c     -----------
435         
436         REAL T
437         DOUBLE PRECISION x
438         DOUBLE PRECISION rho !kg.m-3
439
440c     outputs
441c     -----------
442
443         DOUBLE PRECISION Kthmix
444
445c     local
446c    ------------
447
448         DOUBLE PRECISION x1,x2
449
450         DOUBLE PRECISION  Tc1, Tc2, Pc1, Pc2
451
452         DOUBLE PRECISION  A12, A11, A22, A21
453
454         DOUBLE PRECISION  Gamma1, Gamma2, M1, M2
455         DOUBLE PRECISION  lambda_trans1, lambda_trans2,epsilon
456
457         DOUBLE PRECISION  kco2, kn2
458
459
460      x1 = x
461      x2 = 1d0 - x
462
463
464      M1 = mco2
465      M2 = mn2
466
467
468      Tc1 =  304.1282 !(Scalabrin et al. 2006)
469      Tc2 =  126.192  ! (Lemmon & Jacobsen 2003)
470
471      Pc1 =  73.773   ! (bars)
472      Pc2 =  33.958   ! (bars)
473 
474 
475      Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.)
476      Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.)
477
478
479c Translational conductivities
480
481
482
483      lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) )
484     &                                                          /Gamma1
485
486      lambda_trans2 = ( exp(0.0464 * T/Tc2) - exp(-0.2412 * T/Tc2) )
487     &                                                          /Gamma2
488     
489         
490c     Coefficient of Mason and Saxena
491
492
493      epsilon = 1.
494
495
496      A11 = 1.
497         
498      A22 = 1.
499
500
501      A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)*
502     &                    (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2))
503
504      A21 = epsilon * (1. + sqrt(lambda_trans2/lambda_trans1)*
505     &                    (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1))
506
507
508c     INDIVIDUAL COND.
509
510         call KthCO2Scalab(kco2,T,rho)
511         call KthN2LemJac(kn2,T,rho)
512
513
514c     MIXTURE COND.
515
516        Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22)
517        Kthmix = Kthmix*1e-3   ! from mW.m-1.K-1 to  W.m-1.K-1
518
519        RETURN
520
521        END
522
523
524c======================================================================
525
526         subroutine KthN2LemJac(kthn2,T,rho)
527
528c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
529cWITH viscosity
530c======================================================================
531
532       implicit none
533
534        include "microphys.h"
535c        include "microphysCO2.h"
536
537
538c      arguments
539c     -----------
540         
541         REAL T
542         DOUBLE PRECISION rho !kg.m-3
543
544c     outputs
545c     -----------
546
547         DOUBLE PRECISION kthn2
548
549c     local
550c    ------------
551
552        DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
553        DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
554        DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
555        DOUBLE PRECISION d4,d5,d6,d7,d8,d9
556        DOUBLE PRECISION l4,l5,l6,l7,l8,l9
557        DOUBLE PRECISION t2,t3,t4,t5,t6,t7,t8,t9
558        DOUBLE PRECISION gamma4,gamma5,gamma6,gamma7,gamma8,gamma9
559
560        DOUBLE PRECISION Tc,rhoc
561
562        DOUBLE PRECISION tau, delta
563
564        DOUBLE PRECISION visco
565
566        DOUBLE PRECISION k1, k2
567
568
569         N1 = 1.511d0
570         N2 = 2.117d0
571         N3 = -3.332d0
572
573         N4 = 8.862
574         N5 = 31.11
575         N6 = -73.13
576         N7 = 20.03
577         N8 = -0.7096
578         N9 = 0.2672
579
580         t2 = -1.0d0
581         t3 = -0.7d0
582         t4 = 0.0d0
583         t5 = 0.03
584         t6 = 0.2
585         t7 = 0.8
586         t8 = 0.6
587         t9 = 1.9
588   
589         d4 =  1.
590         d5 =  2.
591         d6 =  3.
592         d7 =  4.
593         d8 =  8.
594         d9 = 10.
595   
596         l4 = 0.
597         gamma4 = 0.
598       
599         l5 = 0.
600         gamma5 = 0.
601       
602         l6 = 1.
603         gamma6 = 1.
604
605         l7 = 2.
606         gamma7 = 1.
607
608         l8 = 2.
609         gamma8 = 1.
610
611         l9 = 2.
612         gamma9 = 1.
613
614
615
616
617c---------------------------------------------------------------------- 
618         
619         call viscoN2(T,visco)  !! v given in microPa.s
620
621       
622         Tc   = 126.192d0
623         rhoc = 11.1839  * 1000 * mn2   !!!from mol.dm-3 to kg.m-3
624
625         tau  = Tc / T
626         delta = rho/rhoc
627
628
629
630         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1
631     
632
633c--------- residual thermal conductivity
634
635
636
637         k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4)         
638     &  +     N5 * tau**t5 * delta**d5 * exp(-gamma5*delta**l5)
639     &  +     N6 * tau**t6 * delta**d6 * exp(-gamma6*delta**l6)
640     &  +     N7 * tau**t7 * delta**d7 * exp(-gamma7*delta**l7)
641     &  +     N8 * tau**t8 * delta**d8 * exp(-gamma8*delta**l8)
642     &  +     N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9)
643
644
645         kthn2 = k1 + k2
646
647
648         RETURN
649
650         END
651
652
653c======================================================================
654
655         subroutine viscoN2(T,visco)
656
657c        Compute viscosity of N2 (Lemmon and Jacobsen, 2003)
658
659c======================================================================
660
661         implicit none
662
663       include "microphys.h"
664c       include "microphysCO2.h"
665c      arguments
666c     -----------
667         
668      REAL T
669
670c     outputs
671c     -----------
672
673      DOUBLE PRECISION visco
674
675
676c     local
677c    ------------
678
679      DOUBLE PRECISION a0,a1,a2,a3,a4
680      DOUBLE PRECISION Tstar,factor,sigma,M2
681      DOUBLE PRECISION RGCS
682     
683
684c---------------------------------------------------------------------- 
685
686 
687      factor = 98.94   ! (K)   
688 
689      sigma  = 0.3656  ! (nm)
690 
691      a0 =  0.431
692      a1 = -0.4623
693      a2 =  0.08406
694      a3 =  0.005341
695      a4 = -0.00331
696 
697      M2 = mn2 * 1.00e3   !!! to g.mol-1
698 
699      Tstar = T*1./factor
700
701      RGCS = exp( a0 + a1 * log(Tstar) + a2 * (log(Tstar))**2. +
702     &                a3 * (log(Tstar))**3. + a4 * (log(Tstar))**4. )
703 
704 
705      visco = 0.0266958 * sqrt(M2*T) / ( sigma**2. * RGCS )  !!! microPa.s
706
707
708      RETURN
709
710      END
711
712
713c======================================================================
714
715         subroutine KthCO2Scalab(kthco2,T,rho)
716
717c        Compute thermal cond of CO2 (Scalabrin et al. 2006)
718
719c======================================================================
720
721         implicit none
722
723
724
725c      arguments
726c     -----------
727
728      REAL T
729      DOUBLE PRECISION rho
730
731c     outputs
732c     -----------
733
734      DOUBLE PRECISION kthco2
735
736c     LOCAL
737c     -----------
738
739      DOUBLE PRECISION Tc,Pc,rhoc, Lambdac
740     
741      DOUBLE PRECISION Tr, rhor, k1, k2
742
743      DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
744      DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
745      DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
746
747     
748
749      Tc   = 304.1282   !(K)
750      Pc   = 7.3773e6   !(MPa)
751      rhoc = 467.6      !(kg.m-3)
752      Lambdac = 4.81384 !(mW.m-1K-1)
753 
754      g1 = 0.
755      g2 = 0.
756      g3 = 1.5
757      g4 = 0.0
758      g5 = 1.0
759      g6 = 1.5
760      g7 = 1.5
761      g8 = 1.5
762      g9 = 3.5
763      g10 = 5.5
764
765
766      h1 = 1.
767      h2 = 5.
768      h3 = 1.
769      h4 = 1.
770      h5 = 2.
771      h6 = 0.
772      h7 = 5.0
773      h8 = 9.0
774      h9 = 0.
775      h10 = 0.
776
777      n1 = 7.69857587
778      n2 = 0.159885811
779      n3 = 1.56918621
780      n4 = -6.73400790
781      n5 = 16.3890156
782      n6 = 3.69415242
783      n7 = 22.3205514
784      n8 = 66.1420950
785      n9 = -0.171779133
786      n10 = 0.00433043347
787
788
789
790      Tr   = T/Tc
791      rhor = rho/rhoc
792
793
794
795      k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2)
796     &     + n3*Tr**(g3)*rhor**(h3)
797
798      k2 = n4*Tr**(g4)*rhor**(h4) + n5*Tr**(g5)*rhor**(h5) 
799     &    + n6*Tr**(g6)*rhor**(h6) + n7*Tr**(g7)*rhor**(h7)
800     &    + n8*Tr**(g8)*rhor**(h8) + n9*Tr**(g9)*rhor**(h9)
801     &    + n10*Tr**(g10)*rhor**(h10)
802   
803      k2  = exp(-5.*rhor**(2.)) * k2
804       
805       
806      kthco2 = (k1 + k2) *  Lambdac   ! mW
807
808
809      RETURN
810
811      END
Note: See TracBrowser for help on using the repository browser.