subroutine evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev, & u,v,tsurf,pqCH4, & tankCH4,pdqCH4,dtsurfevap) use radcommon_h, only: gzlat ! Gravity of the planet [m.s-2] use planete_mod, only: z0 ! Surface roughness of the planet [m] use geometry_mod, only: latitude_deg ! Latitude grid of the planet [°] implicit none !==================================================================== ! ! Purpose ! ------- ! Surface flux for methane evaporation. ! The routine calculates the surface flux of methane ! And the latente heat of methane evaporation. ! ! Implicit scheme. ! We start by adding to variables x the physical tendencies already computed. ! We resolve the equation : ! x(t+1) = x(t) + A * (dx/dt) * dt ! ! ! INPUT ! ----- ! ngrid = Number of grid points in the chunk [-] ! nlayer = Number of vertical layers [-] ! ptimestep = Time step [s] ! pplev = Intermediate pressure levels [Pa] ! zzlay = Altitude at the middle of the atmospheric layers (ref : local surf) [m] ! zzlev = Altitude at the atmospheric layer boundaries (ref : local surf) [m] ! u = Zonal component of the wind [m.s-1] ! v = Meridional component of the wind [m.s-1] ! tsurf = Surface temperature [K] ! pqCH4 = Molar fraction of CH4 [mol.mol-1] ! ! ! OUTPUT ! ------ ! tankCH4 = Depth of surface methane tank [m] ! pdqCH4 = Molar fraction tendency of CH4 at the surface [mol.mol-1.s-1] ! dtsurfevap = Evaporation heating rate at the surface [K.s-1] ! ! ! Author ! ------ ! B. de Batz de Trenquelléon (10/2022) ! !==================================================================== !------------------------------------ ! 0. DECLARATIONS !------------------------------------ ! Inputs : !--------- integer, intent(in) :: ngrid integer, intent(in) :: nlayer real, intent(in) :: ptimestep real, intent(in) :: pplev(ngrid,nlayer+1) real, intent(in) :: zzlay(ngrid,nlayer) real, intent(in) :: zzlev(ngrid,nlayer+1) real, intent(in) :: u(ngrid,nlayer) real, intent(in) :: v(ngrid,nlayer) real, intent(in) :: tsurf(ngrid) real, intent(in) :: pqCH4(ngrid,nlayer) ! Outputs : !---------- real, intent(out) :: tankCH4(ngrid) real, intent(out) :: pdqCH4(ngrid) real, intent(out) :: dtsurfevap(ngrid) ! Parameters : !------------- real, parameter :: karman = 0.4 ! Karman constant [-] real, parameter :: humCH4 = 0.4 ! Imposed surface humidity for CH4 [-] real, parameter :: Flnp = 0.08 ! Fraction occupied by lakes (North Pole) real, parameter :: Flsp = 0.02 ! Fraction occupied by lakes (South Pole) real, parameter :: Flml = 0.02 ! Fraction not infiltrated into the ground (Mid latitudes) real, parameter :: mmolair = 28.e-3 ! Molar mass of air [kg.mol-1] real, parameter :: mmolCH4 = 16.e-3 ! Molar mass of CH4 [kg.mol-1] real, parameter :: rhoiCH4 = 425. ! Density of ice of CH4 [kg.m-3] real, parameter :: TcCH4 = 190.56 ! Critical point of CH4 [K] real, parameter :: Cplake = 2689992 ! Surface heat capacity for hydrocarbon lakes [J.m-2.K-1] (Tokano 2005) ! Local variables : !------------------ integer :: ig ! Initialisation : real*8 :: Tlake ! Lake temperature [K] real*8 :: rhoair ! Density of air [kg.m-3] real*8 :: ws ! Horizontal wind at the surface [m.s-1] real*8 :: Cd ! Turbulent term [-] real*8 :: qsatCH4 ! Saturation profile of CH4 [mol.mol-1] ! Flux : real*8 :: flux ! Surface flux [kg.m-2.s-1] real*8 :: fluxCH4 ! Surface flux fo CH4 [kg.m-2.s-1.mol.mol-1] ! Variations of CH4 : real*8 :: newpqCH4 ! New molar fraction of CH4 in the first layer [mol.mol-1] real*8 :: dtankCH4 ! Trend of CH4's tank [m] ! Latente heat : real*8 :: ftm, LvCH4 ! Variables for latente heat [-, J.kg-1] !------------------------------------ ! 1. INITIALISATION !------------------------------------ do ig = 1, ngrid ! Main loop on the horizontal grid ! Density of the first layer of the atmosphere [kg.m-3] rhoair = (pplev(ig,1) - pplev(ig,2)) / gzlat(ig,1) / (zzlev(ig,2) - zzlev(ig,1)) ! Horizontal winds at the surface [m.s-1] ws = sqrt(u(ig,1)*u(ig,1) + v(ig,1)*v(ig,1)) * (10. / zzlay(ig,1))**0.2 ! Source of turbulent kinetic energy at the surface [-] (Forget et al. 1999) Cd = (karman / log(1. + zzlay(ig,1)/z0))**2 ! Saturation profile of CH4 [mol.mol-1] (Fray and Schmidt 2009) : Tlake = tsurf(ig) - 7 ! Lakes are 2-7K less than surface. qsatCH4 = (1.0e5 / pplev(ig,1)) * exp(1.051e1 - 1.110e3/Tlake - 4.341e3/Tlake**2 + 1.035e5/Tlake**3 - 7.910e5/Tlake**4) ! CH4 : 0.80 * qsat because of dissolution in N2 qsatCH4 = 0.80 * qsatCH4 ! Flux at the surface [kg.m-2.s-1] : flux = rhoair * Cd * ws ! Surface humidity : qsatCH4 = humCH4 * qsatCH4 ! if (REAL(latitude_deg(ig)) .ge. 70.0) then flux = Flnp * flux ! else if (REAL(latitude_deg(ig)) .le. -70.0) then flux = Flsp * flux ! else flux = Flml * flux endif ! Empty tank ? if (tankCH4(ig) .le. 1.e-30) then flux = 0.0 tankCH4(ig) = 1.e-30 endif !------------------------------------ ! 2. IMPLICIT SCHEME !------------------------------------ ! Flux of CH4 at the surface [kg.m-2.s-1.mol.mol-1] : fluxCH4 = flux * (qsatCH4 - pqCH4(ig,1)) ! Flux at the surface [kg.m-2.s-1] --> [s-1] : flux = flux / rhoair / (zzlev(ig,2) - zzlev(ig,1)) ! New molar fraction of CH4 in the first layer [mol.mol-1] : newpqCH4 = (pqCH4(ig,1) + flux * ptimestep * qsatCH4) / (1. + flux * ptimestep) ! Trend of CH4's tank [m] dtankCH4 = - (newpqCH4 - pqCH4(ig,1)) * rhoair * (zzlev(ig,2) - zzlev(ig,1)) * mmolCH4 / mmolair / rhoiCH4 ! New tank depth of CH4 [m] if ((tankCH4(ig) + dtankCH4) .ge. 1.e-30) then tankCH4(ig) = tankCH4(ig) + dtankCH4 else !write(*,*) 'Evaporation CH4 : Empty lakes...' newpqCH4 = tankCH4(ig) / (rhoair * (zzlev(ig,2) - zzlev(ig,1)) * mmolCH4 / mmolair / rhoiCH4) + pqCH4(ig,1) tankCH4(ig) = 1.e-30 endif ! Trend of CH4's molar fraction in the first layer [mol.mol-1.s-1] pdqCH4(ig) = (newpqCH4 - pqCH4(ig,1)) / ptimestep !------------------------------------ ! 3. LATENTE HEAT !------------------------------------ ftm = (1. - Tlake / TcCH4) if(ftm .le. 1.e-3) then ftm = 1.e-3 endif ! Latente heat of CH4 vaporisation [J.kg-1.mol.mol-1] LvCH4 = 8.314 * TcCH4 * (7.08 * ftm**0.354 + 10.95 * 1.1e-2 * ftm**0.456) / mmolCH4 ! Evaporation heating rate [K.s-1] dtsurfevap(ig) = - (fluxCH4 * LvCH4) / Cplake enddo ! End of main loop on the horizontal grid return end subroutine evapCH4