Ignore:
Timestamp:
Nov 2, 2017, 4:40:47 PM (7 years ago)
Author:
jaudouard
Message:

Commit for CO2 clouds microphysics.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F

    r1720 r1816  
    1 
    2 
    31c=======================================================================
    42      subroutine massflowrateCO2(P,T,Sat,Radius,Matm,Ic)
    53c
    6 c     Determination of the mass transfer rate
     4c     Determination of the mass transfer rate for CO2 condensation &
     5c     sublimation
    76c
    87c     newton-raphson method
    9 
    108c     CLASSICAL  (no SF etc.)
    11      
    12 c     AUTOMATIC SETTING OF RANGES FOR NEWTON-RAPHSON FOR THE PAPER
    13 
    14 c     MASS FLUX Ic
    15 
     9c
     10c   inputs: Pressure (P), Temperature (T), saturation ratio (Sat),
     11c           particle radius (Radius), molecular mass of the atm (Matm)
     12c   output:  MASS FLUX Ic
     13c
     14c   Authors: C. Listowski (2014) then J. Audouard (2016-2017)
    1615c=======================================================================
    1716      USE comcstfi_h
    … …  
    2019
    2120      include "microphys.h"
    22 c      include "microphysCO2.h"
    2321
    2422
    … …  
    2624c   arguments: INPUT
    2725c   ----------
    28 
    2926      REAL T,Matm
    3027      REAL*8 SAT
    … …  
    3330c   arguments: OUTPUT
    3431c   ----------
    35 
    3632      DOUBLE PRECISION   Ic
    37 
    3833c   Local Variables
    3934c   ----------
    40 
    4135      DOUBLE PRECISION   Tcm
    4236      DOUBLE PRECISION   T_inf, T_sup, T_dT
    … …  
    4539      DOUBLE PRECISION   rtsafe
    4640      DOUBLE PRECISION   left, fval, dfval
    47 
    4841c    function for newton-raphson iterative method
    4942c    --------------------------
    50 
    5143      EXTERNAL classical
    52 
    5344         
    54       Tcm = 80!dble(T)    ! initialize pourquoi 0 et pas t(i)
    55 
     45      Tcm = 110!dble(T)    ! initialize pourquoi 0 et pas t(i)
    5646      T_inf = 0d0
    5747      T_sup = 800d0
    58 
    59       T_dT = 0.01  ! precision - mettre petit et limiter nb iteration?
     48      T_dT = 0.001  ! precision - mettre petit et limiter nb iteration?
    6049     
    6150666   CONTINUE
    6251
    63 c      print*, 'Radius ', Radius
    64 c      print*, 'SAT = ', Sat
    6552      call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2)
    6653
    67       if (isnan(C0) .eqv. .true.)  C0=0d0
    68 
    69    
     54      if (isnan(C0) .eqv. .true.)  C0=0d0   
    7055c     FIND SURFACE TEMPERATURE (Tc) : iteration sur t
    71 
    7256      cond = 4.*pi*Radius*kmix
    73  
    7457      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 
     58      if (Tcm.LE.0d0 .or. isnan(Tcm)) then ! unsignificant cases where S<<<Seq and Ncores <<1e-10
     59         Tcm = T
    8160      endif
    82 
    83 
    8461c     THEN COMPUTE MASS FLUX Ic from FINAL Tsurface (Tc)
    85 
    8662      Ic = (Tcm-T)
    87 
    88       Ic = cond*Ic/(-Lsub)
     63      Ic = cond*Ic/Lsub
    8964c regarder de combien varie la solution Ic entre Tcm et Tcm+T_dT
    90    
     65 
    9166      RETURN
    9267
    … …  
    9570
    9671c****************************************************************
    97 
    9872      FUNCTION rtsafe(funcd,x1,x2,xacc,Radius,C0,C1,C2)
    99 *
    10073*
    10174*     Newton Raphsen routine (Numerical Recipe)
    … …  
    11487      EXTERNAL funcd
    11588
    116       PARAMETER (MAXIT=50000)   !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection,
     89      PARAMETER (MAXIT=10000)   !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection,
    11790                                !find the root of a function bracketed between x1 and x2. The root, returned as the function value rtsafe,
    11891                                !will be refined until its accuracy is known within !!±xacc. funcd is a user-supplied subroutine which
    … …  
    128101
    129102      if ((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.) ) then
    130          
    131          
    132103         x1=0d0
    133          x2=500d0
    134          
     104         x2=1200d0         
    135105         call funcd(x1,fl,df,C0,C1,C2)
    136106         call funcd(x2,fh,df,C0,C1,C2)
    137          
    138          write(*,*) 'root must be bracketed in rtsafe'
     107c         write(*,*) 'root must be bracketed in rtsafe'
    139108      endif
    140109
    141110
    142111      if (fl.eq.0.) then
    143 
    144112         rtsafe=x1
    145113         return
    146 
    147114      else if (fh.eq.0.) then
    148 
    149115         rtsafe=x2
    150116         return
    151 
    152       else if (fl.lt.0.) then   !Orient the search so that f(xl) < 0.
    153                                                                                                
     117      else if (fl.lt.0.) then   !Orient the search so that f(xl) < 0.                                                                   
    154118         xl=x1
    155119         xh=x2
    156 
    157120      else
    158 
    159121        xh=x1
    160122        xl=x2
    161 
    162123      endif
    163 
    164124      rtsafe = .5*(x1+x2)       !Initialize the guess for root,
    165125      dxold  = abs(x2-x1)       !the stepsize before last,
    166126      dx     = dxold            ! and the last step.
    167 
    168 
    169127      call funcd(rtsafe,f,df,C0,C1,C2)
    170 
    171128      DO 11 j=1,MAXIT           !Loop over allowed iterations.
    172 
    173          
    174          !print*, 'iteration:', j
    175          !print*, rtsafe
    176 
    177129
    178130         if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).gt.0.   ! Bisect if Newton out of range
    … …  
    182134             dx=0.5*(xh-xl)
    183135             rtsafe=xl+dx
    184 
    185136             if (xl.eq.rtsafe) return                   !Change in root is negligible. Newton step acceptable. Take it.
    186 
    187137         else
    188 
    189138            dxold=dx
    190139            dx=f/df
    191140            temp=rtsafe
    192 
    193141            rtsafe=rtsafe-dx
    194 
    195142            if(temp.eq.rtsafe) return
    196 
    197143         endif
    198144
    199 
    200145        if(abs(dx).lt.xacc) return !Convergence criterion. The one new function evaluation per iteration. Maintain the bracket on the root.
    201 
    202146         call funcd(rtsafe,f,df,C0,C1,C2)
    203 
    204147        if(f.lt.0.) then
    205148         xl=rtsafe
    … …  
    207150         xh=rtsafe
    208151        endif
    209 
    210 
    21115211    ENDDO
    212153
    213       write(*,*) 'rtsafe exceeding maximum iterations'
    214 
     154      rtsafe=0d0
     155      write(*,*) 'rtsafe exceeding maximum iterations,Tcm=',rtsafe
    215156      return
    216157
    … …  
    219160
    220161c********************************************************************************
    221 
    222162      subroutine classical(x,f,df,C0,C1,C2)
    223              
    224163c     Function to give as input to RTSAFE (NEWTON-RAPHOEN)
    225 
    226 
    227164c********************************************************************************
    228165
    … …  
    231168      DOUBLE PRECISION   x
    232169      DOUBLE PRECISION  C0,C1,C2
    233 
    234170      DOUBLE PRECISION f
    235171      DOUBLE PRECISION  df
    236      
    237 
    238      
    239172
    240173      f   = x + C0*exp(C1*x) - C2   ! start f
    … …  
    242175 
    243176      return
    244 
    245177      END 
    246178
    … …  
    250182
    251183c********************************************************************************
    252 c defini la fonction eq 6 papier 2014
     184c defini la fonction eq 6 papier 2014  (Listowski et al., 2014)
    253185      use tracer_mod, only: rho_ice_co2
    254186      USE comcstfi_h
    255187
    256188      implicit none
    257 
    258189      include "microphys.h"
    259 c      include "microphysCO2.h"
    260    
    261190
    262191c   arguments: INPUT
    … …  
    270199c   local:
    271200c   ------
    272 
    273 
    274201      DOUBLE PRECISION Cpatm,Cpn2,Cpco2
    275202      DOUBLE PRECISION psat, xinf, pco2
    … …  
    298225
    299226c     Equilibirum pressure over a flat surface
    300 
    301227      psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T))  ! (Pa)
    302 
    303228c     Compute transport coefficient
    304 
    305229      pco2 = psat * dble(S)
    306 
    307230c     Latent heat of sublimation if CO2  co2 (J.kg-1)
    308231c     version Azreg_Ainou (J/kg) :
    309 
    310232      l0=595594.     
    311233      l1=903.111     
    … …  
    313235      l3=0.0528288
    314236      l4=-0.000103183
    315  
    316237      Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 *
    317238     &     dble(T)**3 + l4 * dble(T)**4 ! J/kg
    318 
    319239c     atmospheric density
    320 
    321240      rhoatm = dble(P*Matm)/(rgp*dble(T))   ! g.m-3
    322241      rhoatm = rhoatm * 1.00e-3 !kg.m-3
    323 
    324242      call  KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2
    325243      call  Diffcoeff(P, T, Dv) 
    … …  
    328246
    329247c     ----- FS correction for Diff
    330 
    331248      vthco2  = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1
    332 
    333249      knudsen = 3*Dv / (vthco2 * rc)
    334 
    335250      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
    336      
    337251      Dv      = Dv / (1. + lambda * knudsen)
    338        
    339252c     ----- FS correction for Kth
    340 
    341253      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 
    343254      Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1
    344 
    345255      lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/
    346256     &     (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer
    347 
    348257      knudsen = lpmt / rc
    349 
    350258      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
    351 
    352259      kmix    = kmix /  (1. + lambda * knudsen)
    353 
    354260c     --------------------- ASSIGN coeff values for FUNCTION
    355 
    356261      xinf = dble(S) * psat / dble(P)
    357 
    358262      Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) ))
    359  
    360263      C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/
    361264     &     (rgp*dble(T)))
    362  
    363265      C1 = Lsub*mco2/(rgp*dble(T)**2)
    364 
    365266      C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T))
    366      
    367267      RETURN
    368268      END
    … …  
    370270
    371271c======================================================================
    372 
    373272      subroutine Diffcoeff(P, T, Diff)
    374 
    375273c     Compute diffusion coefficient CO2/N2
    376274c     cited in Ilona's lecture - from Reid et al. 1987
    377275c======================================================================
    378 
    379 
    380276       IMPLICIT NONE
    381277
    382278       include "microphys.h"
    383 c       include "microphysCO2.h"
    384279
    385280c      arguments
    … …  
    425320c======================================================================
    426321
    427 
    428322       implicit none
    429323       
    430324       include "microphys.h"
    431 c       include "microphysCO2.h"
    432 
    433325c      arguments
    434326c     -----------
    … …  
    457349         DOUBLE PRECISION  kco2, kn2
    458350
    459 
    460351      x1 = x
    461352      x2 = 1d0 - x
    462353
    463 
    464354      M1 = mco2
    465355      M2 = mn2
    466356
    467 
    468357      Tc1 =  304.1282 !(Scalabrin et al. 2006)
    469358      Tc2 =  126.192  ! (Lemmon & Jacobsen 2003)
    … …  
    471360      Pc1 =  73.773   ! (bars)
    472361      Pc2 =  33.958   ! (bars)
    473  
    474  
     362   
    475363      Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.)
    476364      Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.)
    477365
    478 
    479366c Translational conductivities
    480 
    481 
    482367
    483368      lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) )
    … …  
    487372     &                                                          /Gamma2
    488373     
    489          
    490374c     Coefficient of Mason and Saxena
    491 
    492 
    493375      epsilon = 1.
    494 
    495376
    496377      A11 = 1.
    … …  
    498379      A22 = 1.
    499380
    500 
    501381      A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)*
    502382     &                    (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2))
    … …  
    505385     &                    (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1))
    506386
    507 
    508387c     INDIVIDUAL COND.
    509388
    … …  
    511390         call KthN2LemJac(kn2,T,rho)
    512391
    513 
    514392c     MIXTURE COND.
    515 
    516393        Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22)
    517394        Kthmix = Kthmix*1e-3   ! from mW.m-1.K-1 to  W.m-1.K-1
    518 
    519395        RETURN
    520396
    521397        END
    522398
    523 
    524 c======================================================================
    525 
     399c======================================================================
    526400         subroutine KthN2LemJac(kthn2,T,rho)
    527 
    528401c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
    529402cWITH viscosity
    … …  
    565438
    566439        DOUBLE PRECISION k1, k2
    567 
    568440
    569441         N1 = 1.511d0
    … …  
    612484         gamma9 = 1.
    613485
    614 
    615 
    616 
    617 c---------------------------------------------------------------------- 
    618          
     486c----------------------------------------------------------------------           
    619487         call viscoN2(T,visco)  !! v given in microPa.s
    620 
    621        
     488     
    622489         Tc   = 126.192d0
    623490         rhoc = 11.1839  * 1000 * mn2   !!!from mol.dm-3 to kg.m-3
    … …  
    626493         delta = rho/rhoc
    627494
    628 
    629 
    630          k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1
    631      
    632 
     495         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1     
    633496c--------- residual thermal conductivity
    634 
    635 
    636497
    637498         k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4)         
    … …  
    642503     &  +     N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9)
    643504
    644 
    645505         kthn2 = k1 + k2
    646 
    647506
    648507         RETURN
    … …  
    744603      DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
    745604      DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
    746 
    747      
    748605
    749606      Tc   = 304.1282   !(K)
    … …  
    763620      g10 = 5.5
    764621
    765 
    766622      h1 = 1.
    767623      h2 = 5.
    … …  
    786642      n10 = 0.00433043347
    787643
    788 
    789 
    790644      Tr   = T/Tc
    791645      rhor = rho/rhoc
    792 
    793 
    794646
    795647      k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2)
    … …  
    802654   
    803655      k2  = exp(-5.*rhor**(2.)) * k2
    804        
    805        
     656               
    806657      kthco2 = (k1 + k2) *  Lambdac   ! mW
    807658
    808 
    809659      RETURN
    810660
Note: See TracChangeset for help on using the changeset viewer.