atmosphere.f90

!>@Master thesis Pierre Boissinot
!>@date 2017
!>update speciation model 2019 - Lena Noack, Claire Guimond, Gianluigi Ortenzi
module atmosphere
  use precision
  use parameters ! add here parameter sets that are used in subroutines and functions below
  use thermal, only:radiogenic_heating
  use mesher ! pas besoin de donner de donner de precisioncomme only ...
  use particles
  use gas_speciation
  use variables  ! attention il faut prendre le nom du module   
 implicit none

contains

  ! 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

  ! 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)
    ! Input values; intent(in): constant input; intent(inout): can be changed here
    real(dp),intent(in)::t,dt,depl_av,T_mean
    type(variables_unknowns),intent(in)::field ! for T
    type(param_F),intent(in)::pF
    type(param_G),intent(in)::pG
    type(param_I),intent(in)::pI
    type(param_M),intent(in)::pM
    type(param_N),intent(in)::pN
    type(param_P),intent(in)::pP
    type(param_X),intent(in)::pX
    type(particle),intent(inout)::part
    type(mesh_cp),intent(inout)::mesh
    type(cell_properties),intent(inout)::cell
    type(melt_properties),intent(inout)::mc
    integer,intent(in)::di
    real(dp),intent(in)::Psurf

    ! Define local variables
    real(dp)::pi__,FD ! FD correspond to Final Dimensionalization
    real(dp):: KI, KII, X_carbonate_melt , fO2
    integer::i,j,k,m
    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

    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)
   
    ! output per time step
    open(unit=60,file="data_atmosphere.res",form='formatted',status='unknown',POSITION='APPEND')

    pi__ = acos(-1.0_dp) 
    FD=(pF%D**3)*pF%rho ! FD correspond to Final Dimensionalization
   
    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 |&
     & 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

     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,&
!       & 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
       & (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
       & mc%frac_H2O,mc%frac_CO2,mc%frac_H2,mc%frac_CO ! mole fractions

    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,&
              &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%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)

  end subroutine



! 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)
      real(dp),intent(in)::t                                                                
      type(variables_unknowns),intent(in)::field ! for T
      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
      type(cell_properties),intent(inout)::cell
      type(melt_properties),intent(inout)::mc ! stored DeltaF and all atmosphere variables (variable.f90/melt_properties (mc))
      real(dp),intent(in)::Psurf_ND
      type(param_F),intent(in)::pF
      type(param_G),intent(in)::pG
      type(param_N),intent(in)::pN
      type(param_M),intent(in)::pM
      type(param_P),intent(in)::pP
      type(param_X),intent(in)::pX
      integer,intent(in)::di

      real(dp)::KI, KII, X_carbonate_melt,T_K,Tm_surf,P_bar,Psurf,fO2,D,pi__
      real(dp)::frac_CO2,frac_CO,frac_H2O,frac_H2,mole_frac_H2O,mole_frac_CO2,vol
      real(dp)::M_H2O,M_CO2,M_H2,M_CO
      integer::i,j,k,m,jmin,jmax,y
      pi__ = acos(-1.0_dp) !3.1415
      Psurf = Psurf_ND*pF%g*pF%rho*pF%D*(10.0_dp**(-5)) ! Surface pressure in bar 

      if (di.eq.2) then  
        jmin = 1
        jmax = 1
        else            
        jmin = 0
        jmax = mesh%ny+1
      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)

          ! H2O equation ! CHECK: better to use aggregate formulation instead of batch melting?
          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))  
              
          ! 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)

          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
          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)
          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
          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 !
      !!!!!!!!!!!!!
      
      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

      mc%mass_H2O=0.0_dp
      mc%mass_CO2=0.0_dp 
      mc%mass_H2=0.0_dp
      mc%mass_CO=0.0_dp

      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

      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

      do i=0,mesh%nl+1
        do j=jmin,jmax 
          do k=0,mesh%nr+1
            if (mc%DeltaF(i,j,k).gt.0.0_dp) then
              ! 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 
                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

                ! 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)
                call spec_COH_Fegley2013(mole_frac_H2O,mole_frac_CO2,pX%fO2,Psurf,Tm_surf,frac_CO2,&
                                        &frac_CO,frac_H2O,frac_H2)
                ! frac_CO2+frac_CO+frac_H2O+frac_H2 is equal to one
                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)

                ! 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
                  mc%mass_H2O = mc%mass_H2O + (cell%X_H2O(i,j,k)/M_H2O * frac_H2O/mole_frac_H2O * M_H2O)&
                                          & * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
                  mc%mass_H2  = mc%mass_H2  + (cell%X_H2O(i,j,k)/M_H2O * frac_H2/mole_frac_H2O * M_H2)&
                                          & * 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
                  mc%mass_CO2 = mc%mass_CO2 + (cell%X_CO2(i,j,k)/M_CO2 * frac_CO2/mole_frac_CO2 * M_CO2)&
                                          & * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
                  mc%mass_CO  = mc%mass_CO  + (cell%X_CO2(i,j,k)/M_CO2 * frac_CO/mole_frac_CO2 * M_CO)&
                                          & * 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))
              end if
            end if 
          enddo
        enddo
      enddo
      mc%mass_H2O = mc%mass_H2O * D ! in nondimensional Kg
      mc%mass_CO2 = mc%mass_CO2 * D
      mc%mass_H2  = mc%mass_H2  * D
      mc%mass_CO  = mc%mass_CO  * D

      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
       
      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

  end subroutine atmos_outgassing




  subroutine atmos_escaping (t,mc,pF,pM,mesh)
       real(dp),intent(in)::t
       type(param_F),intent(in)::pF
       type(param_M),intent(in)::pM
       type(mesh_cp),intent(inout)::mesh
       type(melt_properties),intent(inout)::mc

       integer::i,j,k,jmin,jmax
       real(dp)::FD,time_escape
!!!!!! initial parameters

   
  !Na= 6.02_dp**23
  !MM_H2O = 18.01_dp
  !MM_CO2 = 44.01_dp
   mc%H2O_atm=0.0_dp
   mc%CO2_atm=0.0_dp
   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

!!!!!! 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 
              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 
                                       ! (part_H2O_atm= (H2O_atm+initial_H2O)/m_part_H2O)
               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 
                                       ! (part_CO2_atm= (CO2_atm+ initial_CO2)/m_part_CO2)
   

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

              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
               
!!!!!! Calculation of the variation of the Atmosphere composition 

              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
               
             !  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

              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
                end if
               end if
             enddo
            enddo
          enddo

              mc%tot_esc_H2O= mc%tot_esc_H2O + mc%escape_H2O
              mc%tot_esc_CO2= mc%tot_esc_CO2 + mc%escape_CO2

 end subroutine atmos_escaping


end module atmosphere