atmosphere.f90 21.2 KB
Newer Older
Lena HPC's avatar
Lena HPC committed
1
2
!>@Master thesis Pierre Boissinot
!>@date 2017
3
!>update speciation model 2019 - Lena Noack, Claire Guimond, Gianluigi Ortenzi
Lena HPC's avatar
Lena HPC committed
4
5
module atmosphere
  use precision
Boissinot Pierre's avatar
Boissinot Pierre committed
6
  use parameters ! add here parameter sets that are used in subroutines and functions below
Lena HPC's avatar
Lena HPC committed
7
  use thermal, only:radiogenic_heating
Boissinot Pierre's avatar
Boissinot Pierre committed
8
9
  use mesher ! pas besoin de donner de donner de precisioncomme only ...
  use particles
10
  use gas_speciation
Boissinot Pierre's avatar
Boissinot Pierre committed
11
12
  use variables  ! attention il faut prendre le nom du module   
 implicit none
Lena HPC's avatar
Lena HPC committed
13
14
15

contains

16
17
18
19
  ! ToDo: needs to be input: pX = param_X (all concentration and atmosphere parameters)
  ! includes: pX%fO2, pX%CO2_ini, pX%fwm, pX%use_atmos (need to be set in input.txt)
  ! available now in part/cell: part%CO2(m), part%X_CO2(m), cell%CO2(i,j,k), cell%X_CO2(i,j,k) -> solid and melt concentrations

20
21
  ! all time-dependent variables like mass_H2O etc. are stored in the mc structure, defined in variables.f90
  subroutine get_atmos(t,dt,pF,pG,pI,pM,pN,pP,pX,field,depl_av,T_mean,part,mesh,cell,mc,di,Psurf)
Lena HPC's avatar
Lena HPC committed
22
23
    ! Input values; intent(in): constant input; intent(inout): can be changed here
    real(dp),intent(in)::t,dt,depl_av,T_mean
Boissinot Pierre's avatar
Boissinot Pierre committed
24
    type(variables_unknowns),intent(in)::field ! for T
Lena HPC's avatar
Lena HPC committed
25
26
27
28
    type(param_F),intent(in)::pF
    type(param_G),intent(in)::pG
    type(param_I),intent(in)::pI
    type(param_M),intent(in)::pM
Boissinot Pierre's avatar
Boissinot Pierre committed
29
30
    type(param_N),intent(in)::pN
    type(param_P),intent(in)::pP
31
    type(param_X),intent(in)::pX
Boissinot Pierre's avatar
Boissinot Pierre committed
32
    type(particle),intent(inout)::part
Boissinot Pierre's avatar
Boissinot Pierre committed
33
    type(mesh_cp),intent(inout)::mesh
Boissinot Pierre's avatar
Boissinot Pierre committed
34
    type(cell_properties),intent(inout)::cell
Boissinot Pierre's avatar
Boissinot Pierre committed
35
    type(melt_properties),intent(inout)::mc
36
37
    integer,intent(in)::di
    real(dp),intent(in)::Psurf
38

39
    ! Define local variables
Boissinot Pierre's avatar
Boissinot Pierre committed
40
    real(dp)::pi__,FD ! FD correspond to Final Dimensionalization
41
42
    real(dp):: KI, KII, X_carbonate_melt , fO2
    integer::i,j,k,m
Claire Guimond's avatar
Claire Guimond committed
43
    real(dp)::M_tot,X_mol_H2,X_mol_H2O,X_mol_CO,X_mol_CO2,sum_moles,frac_tot_H2,frac_tot_H2O,frac_tot_CO,frac_tot_CO2
44
45
46

    call atmos_outgassing(t,part,mesh,cell,mc,field,di,pF,pG,pM,pN,pP,pX,Psurf)
    if (pX%use_atmos.ge.2) call atmos_escaping(t,mc,pF,pM,mesh)
47
   
Lena HPC's avatar
Lena HPC committed
48
49
50
    ! output per time step
    open(unit=60,file="data_atmosphere.res",form='formatted',status='unknown',POSITION='APPEND')

51
    pi__ = acos(-1.0_dp) 
Boissinot Pierre's avatar
Boissinot Pierre committed
52
53
    FD=(pF%D**3)*pF%rho ! FD correspond to Final Dimensionalization
   
54
55
56
57
    if (pX%gas_spec.ge.1) then
     ! no escape included, but different species outgassed
     if (t.eq.dt) write(60,*) '  Time (Gyr) | T_mantle(K) | maxF | mass H2O |&
     & total H2O | mass CO2 | total CO2 | mass H2 | total H2 | mass CO | total CO |&
Claire Guimond's avatar
Claire Guimond committed
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
     & PH2O(bar) | PCO2(bar) | PH2(bar) | PCO(bar) | EGL_H2O | frac_H2O | frac_CO2 | frac_H2 | frac_CO'

     M_tot = mc%total_mass_H2O + mc%total_mass_CO2 + mc%total_mass_H2 + mc%total_mass_CO
     X_mol_H2 = mc%total_mass_H2*1000.0_dp / 2.01588_dp           ! number of moles
     X_mol_H2O = mc%total_mass_H2O*1000.0_dp / 18.01488_dp
     X_mol_CO = mc%total_mass_CO*1000.0_dp / 28.0097_dp
     X_mol_CO2 = mc%total_mass_CO2*1000.0_dp / 44.0087_dp

     sum_moles = X_mol_H2+X_mol_H2O+X_mol_CO+X_mol_CO2 

     if (sum_moles.gt.0.0_dp) then
       frac_tot_H2 = X_mol_H2 / sum_moles
       frac_tot_H2O = X_mol_H2O / sum_moles
       frac_tot_CO = X_mol_CO / sum_moles
       frac_tot_CO2 = X_mol_CO2 / sum_moles
     else
       frac_tot_H2 = 0.0_dp
       frac_tot_H2O = 0.0_dp
       frac_tot_CO = 0.0_dp
       frac_tot_CO2 = 0.0_dp
     endif
79
80
81
82
83
84
85
86
87
88
89
90
91
92

     write(60,'(F13.8,F12.3,F13.8,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,&
       & F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6)') &
       & t/pF%time_yr*1.0e-9_dp,& ! time in Gyr
       & T_mean*pF%DeltaT+pF%Ts,& ! average mantle temperature in K
       & maxval(part%df),& ! to see when melting starts
       & mc%mass_H2O*FD,&
       & mc%total_mass_H2O*FD,&
       & mc%mass_CO2*FD,&
       & mc%total_mass_CO2*FD,&
       & mc%mass_H2*FD,&
       & mc%total_mass_H2*FD,&
       & mc%mass_CO*FD,&
       & mc%total_mass_CO*FD,&
Claire Guimond's avatar
Claire Guimond committed
93
94
95
96
97
98
99
100
!       & mc%total_mass_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2O accumulating over time in bar
!       & mc%total_mass_CO2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO2 accumulating over time in bar
!       & mc%total_mass_H2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2 accumulating over time in bar
!       & mc%total_mass_CO*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO accumulating over time in bar
       & M_tot * frac_tot_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2O accumulating over time in bar
       & M_tot * frac_tot_CO2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO2 accumulating over time in bar
       & M_tot * frac_tot_H2 *pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2 accumulating over time in bar
       & M_tot * frac_tot_CO *pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO accumulating over time in bar
101
       & (mc%total_mass_H2O*FD - mc%escape_H2O)/(4.0_dp*pi__*mesh%rmax**2*pF%D**2*1000.0_dp),& ! EGL H2O in m, assuming rho_water=1000  !X_atm_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*1e5_dp/(pF%g*1000.0_dp) ! EGL H2O
Claire Guimond's avatar
Claire Guimond committed
102
       & mc%frac_H2O,mc%frac_CO2,mc%frac_H2,mc%frac_CO ! mole fractions
103
104
105
106
107
108
109
110
111

    else

     if (t.eq.dt) write(60,*) '  Time (Gyr)    T_mantle(K)   maxF          massH2O    &
     &     totalH2O        massCO2    total CO2      PH2O(bar)   PCO2(bar)    max_CO2_frac  &
     & esc_H2O_(kg/Gyr) esc_CO2_(kg/Gyr) tot_H2O_atm     tot_CO2_atm  PH2O_atm_end PCO2_atm_end &
     &EGL_H2O '

     write(60,'(F13.8,F12.3,F13.8,ES16.4,ES16.4,ES16.4,ES13.6,F12.6,F12.6,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,&
Boissinot Pierre's avatar
Boissinot Pierre committed
112
              &F12.6,F12.6,F12.6)') &
Lena HPC's avatar
Lena HPC committed
113
       & t/pF%time_yr*1.0e-9_dp,& ! time in Gyr
114
       & T_mean*pF%DeltaT+pF%Ts,& ! average mantle temperature in K
115
       & maxval(part%df),& ! to see when melting starts
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
       & mc%mass_H2O*FD,&
       & mc%total_mass_H2O*FD,&
       & mc%mass_CO2*FD,&
       & mc%total_mass_CO2*FD,&
       & mc%total_mass_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2O accumulating over time in bar
       & mc%total_mass_CO2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO2 accumulating over time in bar
       & mc%max_CO2_frac,&
       & mc%escape_H2O,& ! in kg/Gyr
       & mc%escape_CO2,& ! in kg/Gyr
       & mc%total_mass_H2O*FD - mc%escape_H2O,& ! in kg !X_atm_H2O,& ! in kg
       & mc%total_mass_CO2*FD - mc%escape_CO2,& ! in kg
       & (mc%total_mass_H2O*FD - mc%escape_H2O)*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5),& !X_atm_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5),& ! atm pressure with outgassing and escaping parts for H2O
       & (mc%total_mass_CO2*FD - mc%escape_CO2)*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5),& ! atm pressure with outgassing and escaping parts for CO2
       & (mc%total_mass_H2O*FD - mc%escape_H2O)/(4.0_dp*pi__*mesh%rmax**2*pF%D**2*1000.0_dp) ! EGL H2O in m, assuming rho_water=1000  !X_atm_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*1e5_dp/(pF%g*1000.0_dp) ! EGL H2O

    endif
         
    close(60)
Lena HPC's avatar
Lena HPC committed
134
135
136

  end subroutine

Lena HPC's avatar
Lena HPC committed
137

Boissinot Pierre's avatar
Boissinot Pierre committed
138

139
140
! calculate the amount of volatiles being outgassed per time step
  subroutine atmos_outgassing(t,part,mesh,cell,mc,field,di,pF,pG,pM,pN,pP,pX,Psurf_ND)
Boissinot Pierre's avatar
Boissinot Pierre committed
141
      real(dp),intent(in)::t                                                                
Boissinot Pierre's avatar
Boissinot Pierre committed
142
      type(variables_unknowns),intent(in)::field ! for T
Boissinot Pierre's avatar
Boissinot Pierre committed
143
144
      type(particle),intent(inout)::part ! df is like DeltaF but it's define like df in the file particle.
      type(mesh_cp),intent(inout)::mesh
Boissinot Pierre's avatar
Boissinot Pierre committed
145
      type(cell_properties),intent(inout)::cell
146
147
      type(melt_properties),intent(inout)::mc ! stored DeltaF and all atmosphere variables (variable.f90/melt_properties (mc))
      real(dp),intent(in)::Psurf_ND
148
      type(param_F),intent(in)::pF
Boissinot Pierre's avatar
Boissinot Pierre committed
149
150
      type(param_G),intent(in)::pG
      type(param_N),intent(in)::pN
151
      type(param_M),intent(in)::pM
Boissinot Pierre's avatar
Boissinot Pierre committed
152
153
      type(param_P),intent(in)::pP
      type(param_X),intent(in)::pX
154
155
      integer,intent(in)::di

156
      real(dp)::KI, KII, X_carbonate_melt,T_K,Tm_surf,P_bar,Psurf,fO2,D,pi__
Claire Guimond's avatar
Claire Guimond committed
157
      real(dp)::frac_CO2,frac_CO,frac_H2O,frac_H2,mole_frac_H2O,mole_frac_CO2,vol
158
      real(dp)::M_H2O,M_CO2,M_H2,M_CO
Boissinot Pierre's avatar
Boissinot Pierre committed
159
      integer::i,j,k,m,jmin,jmax,y
Boissinot Pierre's avatar
Boissinot Pierre committed
160
      pi__ = acos(-1.0_dp) !3.1415
161
      Psurf = Psurf_ND*pF%g*pF%rho*pF%D*(10.0_dp**(-5)) ! Surface pressure in bar 
Boissinot Pierre's avatar
Boissinot Pierre committed
162

Lena Noack's avatar
Lena Noack committed
163
      if (di.eq.2) then  
Boissinot Pierre's avatar
Boissinot Pierre committed
164
165
        jmin = 1
        jmax = 1
Boissinot Pierre's avatar
Boissinot Pierre committed
166
        else            
Boissinot Pierre's avatar
Boissinot Pierre committed
167
168
        jmin = 0
        jmax = mesh%ny+1
Lena Noack's avatar
Lena Noack committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
      endif

      !!!!!!!!!!!!!!!!!
      ! for particles !
      !!!!!!!!!!!!!!!!!

      ! equation X_H2O 18 from Katz et al., (2003) / equation KI and KII from Holloway et al., (1992) 

      part%X_CO2=0.0_dp
      part%X_H2O=0.0_dp 
      do m=1,part%n 
        if (part%df(m).gt.0.0_dp) then
          i = part%cell_l(m)
          j = part%cell_y(m)
          k = part%cell_r(m)

185
          ! H2O equation ! CHECK: better to use aggregate formulation instead of batch melting?
Lena Noack's avatar
Lena Noack committed
186
187
188
189
190
          part%X_H2O(m) = (part%w(m))/(pX%dH2O+(part%df(m)*(1.0_dp-pX%dH2O))) !-> instead: cell%X_H2O / part%X_H2O
          if (part%X_H2O(m)*part%df(m).gt.part%w(m)) then ! condition for having more water after than before the calculation
            part%X_H2O(m)= part%w(m)/part%df(m) 
          end if
          part%w(m) = part%w(m) - (part%X_H2O(m)*part%df(m))  
Boissinot Pierre's avatar
Boissinot Pierre committed
191
              
Lena Noack's avatar
Lena Noack committed
192
193
194
195
196
197
198
199
200
          ! CO2 equations (Holloway)
          T_K= (field%T(i,j,k)*pF%DeltaT)+pF%Ts
          P_bar= (cell%ref_p(i,j,k)*pF%g*pF%rho*pF%D)*(10.0_dp**(-5)) ! Pressure in bar

          KI= 40.07639_dp-((2.53932_dp*10.0_dp**(-2))*T_K)+(5.27096_dp*10.0_dp**(-6))*(T_K**2)+0.0267_dp*(P_bar-1.0_dp)/T_K
          KII= -6.24763_dp-(282.56_dp/T_K)-(0.119242_dp*(P_bar-1000.0_dp)/T_K)
          KI=10.0_dp**(KI)  ! for the inverse of LOG10
          KII= 10.0_dp**(KII)

201
          fO2= 6.899_dp-(27714.0_dp/T_K)+0.05_dp*(P_bar-1.0_dp)/T_K + pX%fO2 ! instead fO2-field, where with change in O atoms we calculate with the same formula back how fO2 changes locally
Lena Noack's avatar
Lena Noack committed
202
203
204
          fO2= 10.0_dp**(fO2)
          X_carbonate_melt = (KI*KII*fO2)/(1.0_dp+(KI*KII*fO2))  
          part%X_CO2(m)=((44.01_dp/pX%fwm)*X_carbonate_melt)/(1.0_dp+(1.0_dp-(44.01_dp/pX%fwm))* X_carbonate_melt) ! with 44.01 represent the MM of CO2 (12+2*16)
205
          part%X_CO2(m) = (part%X_CO2(m))*(1.0_dp-(1.0_dp-part%df(m))**(1.0_dp/pX%dCO2))/part%df(m) ! fractional melting
Lena Noack's avatar
Lena Noack committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
          part%X_CO2(m) = part%X_CO2(m) * 10.0_dp**6
          if (part%X_CO2(m)*part%df(m).gt.part%CO2(m)) then
            part%X_CO2(m) = part%CO2(m)/part%df(m)
          end if
          part%CO2(m) = part%CO2(m)-(part%X_CO2(m)*part%df(m)) 
               
        end if
      end do
      call part_part2cell(part,mc,mesh,cell,pG%geom,pN%debug,pG%periodic,pN%average,pN%p_average,pN%p_average2,&
           &pN%pr_time,pP%moveRef)
                           

      !!!!!!!!!!!!!
      ! for cells !
      !!!!!!!!!!!!!
Boissinot Pierre's avatar
Boissinot Pierre committed
221
      
Lena Noack's avatar
Lena Noack committed
222
223
      D=((4*pi__)/3)*((mesh%rmax**(3)-mesh%rmin**(3))/((mesh%rmax**(2)-mesh%rmin**(2)))) !D is for the shift from 2D to 3D calculations

224
225
226
227
228
      mc%mass_H2O=0.0_dp
      mc%mass_CO2=0.0_dp 
      mc%mass_H2=0.0_dp
      mc%mass_CO=0.0_dp

Claire Guimond's avatar
Claire Guimond committed
229
230
231
232
233
234
      mc%frac_CO = 0.0_dp
      mc%frac_CO2 = 0.0_dp
      mc%frac_H2 = 0.0_dp
      mc%frac_H2O = 0.0_dp
      vol = 0.0_dp

235
236
237
238
239
240
241
242
243
      M_H2O = 2.0_dp*1.00794_dp + 15.999_dp
      M_CO2 = 12.0107_dp + 2.0_dp*15.999_dp
      M_H2 = 2.0_dp*1.00794_dp
      M_CO = 12.0107_dp + 15.999_dp

      ! trace how much CO2 goes into melt, not used otherwise
      mc%Concentration_CO2_melt = 0.0_dp
      mc%Concentration_CO2_frac = 0.0_dp
      mc%max_CO2_frac= 0.0_dp
Lena Noack's avatar
Lena Noack committed
244
245
246
247

      do i=0,mesh%nl+1
        do j=jmin,jmax 
          do k=0,mesh%nr+1
248
            if (mc%DeltaF(i,j,k).gt.0.0_dp) then
249
250
251
              ! X_H2O etc are in ppm, mass_H2O etc are weigth amounts (corresponding to kg in dimensional version)
              if (pX%gas_spec.ge.1) then
                ! calc adiabatic melt temp reduction due to ascent to surface 
252
253
254
255
                T_K= (field%T(i,j,k)*pF%DeltaT)+pF%Ts
                P_bar= (cell%ref_p(i,j,k)*pF%g*pF%rho*pF%D)*(10.0_dp**(-5)) ! Pressure in bar
                Tm_surf = T_K*exp(-pX%alpha_melt*pF%alpha/(pX%density_melt*pF%rho*pX%Cp_melt*pF%Cp)*(P_bar-Psurf)*10.0**5) 
                !write(*,*) i,j,k,"T=",T_K,",Tm_surf=",Tm_surf,", P=",P_bar
256
257
258
259
260

                ! go from wt-ppm of H2O/CO2 in melt to mole fractions (ppm does not matter here since only ratio counts)
                ! sum of mole_frac_H2O and mole_frac_CO2 needs to be one
                mole_frac_H2O = cell%X_H2O(i,j,k)/M_H2O/(cell%X_H2O(i,j,k)/M_H2O+cell%X_CO2(i,j,k)/M_CO2)
                mole_frac_CO2 = cell%X_CO2(i,j,k)/M_CO2/(cell%X_H2O(i,j,k)/M_H2O+cell%X_CO2(i,j,k)/M_CO2)
Claire Guimond's avatar
Claire Guimond committed
261
262
                call spec_COH_Fegley2013(mole_frac_H2O,mole_frac_CO2,pX%fO2,Psurf,Tm_surf,frac_CO2,&
                                        &frac_CO,frac_H2O,frac_H2)
263
                ! frac_CO2+frac_CO+frac_H2O+frac_H2 is equal to one
Claire Guimond's avatar
Claire Guimond committed
264
265
266
267
268
                mc%frac_CO = mc%frac_CO + frac_CO * mesh%dVi(i,j,k)
                mc%frac_CO2 = mc%frac_CO2 + frac_CO2 * mesh%dVi(i,j,k)
                mc%frac_H2 = mc%frac_H2 + frac_H2 * mesh%dVi(i,j,k)
                mc%frac_H2O = mc%frac_H2O + frac_H2O * mesh%dVi(i,j,k)
                vol = vol + mesh%dVi(i,j,k)
269
270
271
272

                ! calc mass H2O/CO2/H2/CO each in kg (when dimensionalized, here still dimensionless value)
                ! first calculate mole of H2O and CO2, use frac_X to get moles of each species, multiply with molar masses to get final gas masses
                if (mole_frac_H2O.gt.0.0_dp) then
Claire Guimond's avatar
Claire Guimond committed
273
                  mc%mass_H2O = mc%mass_H2O + (cell%X_H2O(i,j,k)/M_H2O * frac_H2O/mole_frac_H2O * M_H2O)&
274
                                          & * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
Claire Guimond's avatar
Claire Guimond committed
275
                  mc%mass_H2  = mc%mass_H2  + (cell%X_H2O(i,j,k)/M_H2O * frac_H2/mole_frac_H2O * M_H2)&
276
277
278
                                          & * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
                end if
                if (mole_frac_CO2.gt.0.0_dp) then
Claire Guimond's avatar
Claire Guimond committed
279
                  mc%mass_CO2 = mc%mass_CO2 + (cell%X_CO2(i,j,k)/M_CO2 * frac_CO2/mole_frac_CO2 * M_CO2)&
280
                                          & * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
Claire Guimond's avatar
Claire Guimond committed
281
                  mc%mass_CO  = mc%mass_CO  + (cell%X_CO2(i,j,k)/M_CO2 * frac_CO/mole_frac_CO2 * M_CO)&
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
                                          & * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
                end if
              else
                ! calc mass H2O/CO2 each in kg without speciation model
                mc%mass_H2O = mc%mass_H2O+cell%X_H2O(i,j,k)*mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)&
                         &*cell%ref_r(i,j,k)*(1.0e-6_dp) ! partial weight 
                mc%mass_CO2 = mc%mass_CO2+cell%X_CO2(i,j,k)*mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)&
                         &*cell%ref_r(i,j,k)*(1.0e-6_dp) ! partial weight 
              end if

              ! trace CO2 in melt
              if (cell%X_CO2(i,j,k).gt.mc%max_CO2_frac) then
                mc%max_CO2_frac = cell%X_CO2(i,j,k)
                mc%Concentration_CO2_melt = mc%Concentration_CO2_melt + (cell%X_CO2(i,j,k)*mc%DeltaF(i,j,k) * mesh%dVi(i,j,k))    
                mc%Concentration_CO2_frac = mc%Concentration_CO2_frac + ((cell%X_CO2(i,j,k)) * mesh%dVi(i,j,k))
Lena Noack's avatar
Lena Noack committed
297
298
              end if
            end if 
Boissinot Pierre's avatar
Boissinot Pierre committed
299
          enddo
Lena Noack's avatar
Lena Noack committed
300
301
        enddo
      enddo
Claire Guimond's avatar
Claire Guimond committed
302
      mc%mass_H2O = mc%mass_H2O * D ! in nondimensional Kg
303
304
305
      mc%mass_CO2 = mc%mass_CO2 * D
      mc%mass_H2  = mc%mass_H2  * D
      mc%mass_CO  = mc%mass_CO  * D
Claire Guimond's avatar
Claire Guimond committed
306
307
308
309
310

      mc%frac_CO = mc%frac_CO / vol
      mc%frac_CO2 = mc%frac_CO2 / vol
      mc%frac_H2 = mc%frac_H2 / vol
      mc%frac_H2O = mc%frac_H2O / vol
311
       
312
313
314
315
316
      mc%total_mass_H2O = mc%total_mass_H2O + mc%mass_H2O*pM%Chi_extr 
      mc%total_mass_CO2 = mc%total_mass_CO2 + mc%mass_CO2*pM%Chi_extr
      mc%total_mass_H2  = mc%total_mass_H2  + mc%mass_H2 *pM%Chi_extr
      mc%total_mass_CO  = mc%total_mass_CO  + mc%mass_CO *pM%Chi_extr

Lena Noack's avatar
Lena Noack committed
317
  end subroutine atmos_outgassing
318
319
320



Boissinot Pierre's avatar
Boissinot Pierre committed
321

322
  subroutine atmos_escaping (t,mc,pF,pM,mesh)
Boissinot Pierre's avatar
Boissinot Pierre committed
323
       real(dp),intent(in)::t
324
325
326
       type(param_F),intent(in)::pF
       type(param_M),intent(in)::pM
       type(mesh_cp),intent(inout)::mesh
Boissinot Pierre's avatar
Boissinot Pierre committed
327
       type(melt_properties),intent(inout)::mc
328
329

       integer::i,j,k,jmin,jmax
330
       real(dp)::FD,time_escape
331
332
333
334
335
336
!!!!!! initial parameters

   
  !Na= 6.02_dp**23
  !MM_H2O = 18.01_dp
  !MM_CO2 = 44.01_dp
337
338
   mc%H2O_atm=0.0_dp
   mc%CO2_atm=0.0_dp
Boissinot Pierre's avatar
Boissinot Pierre committed
339
340
341
342
   time_escape = t/pF%time_yr*1.0e-9_dp ! variation of time in Gyr with pF%time correspond to the transition from second to year.
   !tot_esc_H2O=0.0_dp
   !tot_esc_CO2=0.0_dp
   FD=(pF%D**3)*pF%rho ! FD correspond to Final Dimensionalization
343
344
345
346

!!!!!! calculation of the real amount of H2O and CO2 in the atm and the particles amount equivalent
         do i=0,mesh%nl+1
            do j=jmin,jmax 
347
348
349
350
351
352
              do k=0,mesh%nr+1         ! ERROR? Chi_extr already used above... Check everywhere in escape routine
               mc%H2O_atm = mc%H2O_atm + (mc%total_mass_H2O*pM%Chi_extr) ! amount of H2O in the atm which is taking into account the extrusive H2O factor following Grott et al., 2011
               mc%CO2_atm = mc%CO2_atm + (mc%total_mass_CO2*pM%Chi_extr) ! amount of CO2 in the atm which is taking into account the extrusive CO2 factor following Grott et al., 2011
               mc%mass_part_H2O = (1.0_dp/(6.02e23_dp)*18.01_dp)*1e-3_dp ! calculation of the mass of one particle of H2O ((1/Na)*MM) in kg and *10.0_dp**(-3) for kg  ((1.0_dp/6.02_dp**(23))*18.01_dp)
               mc%mass_part_CO2 = (1.0_dp/(6.02e23_dp)*44.01_dp)*1e-3_dp  ! calculation of the mass of one particle of CO2
               mc%part_H2O_atm = mc%H2O_atm/mc%mass_part_H2O ! number of particles of H2O in the atm from outgassing simulation / we can add any primary composition here after 
353
                                       ! (part_H2O_atm= (H2O_atm+initial_H2O)/m_part_H2O)
354
               mc%part_CO2_atm = mc%CO2_atm/mc%mass_part_CO2 ! number of particles of CO2 in the atm from outgassing simulation / we can add any primary composition here after 
355
356
357
358
359
                                       ! (part_CO2_atm= (CO2_atm+ initial_CO2)/m_part_CO2)
   

!!!!!! initial escape rates from the biblio and calculation of the variation during time

360
361
362
363
364
365
366
367
368
369
370
              mc%max_modern_rate_H2O= (3.82e25_dp)!*3.1536e16_dp)!*mass_part_H2O ! escape rate at max solar intensity from biblio research in kg_H2O*Gyr-1 (it can be for a sum of non-thermal mechanisms or from only one)  
              mc%min_modern_rate_H2O= (1.42e25_dp)!*3.1536e16_dp)!*mass_part_H2O ! escape rate at min solar intensity from biblio research in kg_CO2*Gyr-1 (it can be for a sum of non-thermal mechanisms or from only one)
              mc%max_modern_rate_CO2= (1.9e24_dp)!*3.1536e16_dp)!*mass_part_CO2 ! escape rate at max solar intensity from biblio research in kg_H2O*Gyr-1 (it can be for a sum of non-thermal mechanisms or from only one) 
              mc%min_modern_rate_CO2= (3.9e22_dp)!*3.1536e16_dp)!*mass_part_CO2 ! escape rate at min solar intensity from biblio research in kg_CO2*Gyr-1 (it can be for a sum of non-thermal mechanisms or from only one) 

              mc%escape_variation_H2O= (mc%min_modern_rate_H2O*(EXP(LOG(mc%max_modern_rate_H2O/mc%min_modern_rate_H2O))&
                                       &/2.5_dp*4.5_dp))&
                                    &*(EXP((-LOG(mc%max_modern_rate_H2O/mc%min_modern_rate_H2O))/2.5_dp*time_escape))! equation of the variation of the escape rate which is calculate from the modern escape rates for max an min solar intensity and in kg
              mc%escape_variation_CO2= (mc%min_modern_rate_CO2*(EXP(LOG(mc%max_modern_rate_CO2/mc%min_modern_rate_CO2))&
                                       &/2.5_dp*4.5_dp))&
                                    &*(EXP((-LOG(mc%max_modern_rate_CO2/mc%min_modern_rate_CO2))/2.5_dp*time_escape))! equation of the variation of the escape rate which is calculate from the modern escape rates for max an min solar intensity and in kg
Boissinot Pierre's avatar
Boissinot Pierre committed
371
               
372
373
!!!!!! Calculation of the variation of the Atmosphere composition 

374
375
              mc%escape_H2O = mc%escape_variation_H2O*(2.9916e-26_dp*1.0e9_dp) ! in kg/Gyr
              mc%escape_CO2 = mc%escape_variation_CO2*(7.19117e-26_dp*1.0e9_dp) ! in kg/Gyr
Boissinot Pierre's avatar
Boissinot Pierre committed
376
377
378
379
380
381
               
             !  if (mc%DeltaF(i,j,k).gt.0.0_dp) then
             !   tot_esc_H2O= tot_esc_H2O + escape_H2O
             !   tot_esc_CO2= tot_esc_CO2 + escape_CO2
             !  end if

382
383
384
385
386
387
              mc%X_atm_H2O= (mc%total_mass_H2O*FD*0.4)- (mc%escape_H2O) !*(t/pF%time_yr*1.0e-9_dp)) ! variation of the mass of H2O in the atm without taking into account any primary composition 
               if(mc%X_atm_H2O.le.0.0_dp) then
                mc%X_atm_H2O=0.0_dp
              mc%X_atm_CO2= (mc%total_mass_CO2*FD*0.4)- (mc%escape_CO2) !*(t/pF%time_yr*1.0e-9_dp)) ! variation of the mass of CO2 in the atm without taking into account any primary composition 
               if(mc%X_atm_CO2.le.0.0_dp) then
                 mc%X_atm_CO2=0.0_dp
Boissinot Pierre's avatar
Boissinot Pierre committed
388
389
                end if
               end if
390
391
392
393
             enddo
            enddo
          enddo

394
395
              mc%tot_esc_H2O= mc%tot_esc_H2O + mc%escape_H2O
              mc%tot_esc_CO2= mc%tot_esc_CO2 + mc%escape_CO2
Boissinot Pierre's avatar
Boissinot Pierre committed
396

397
 end subroutine atmos_escaping
Boissinot Pierre's avatar
Boissinot Pierre committed
398
399


Lena HPC's avatar
Lena HPC committed
400
end module atmosphere
Lena HPC's avatar
Lena HPC committed
401