source: trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90 @ 3083

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(s)
!   ---------
!   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.5     ! Imposed surface humidity for CH4 [-]

real, parameter :: Flnp = 0.01      ! Fraction occupied by lakes (North Pole)
real, parameter :: Flsp = 0.005     ! Fraction occupied by lakes (South Pole)

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  :: 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) :
  qsatCH4 = (1.0e5 / pplev(ig,1)) * exp(1.051e1 - 1.110e3/tsurf(ig) - 4.341e3/tsurf(ig)**2 + 1.035e5/tsurf(ig)**3 - 7.910e5/tsurf(ig)**4)
  qsatCH4 = humCH4 * qsatCH4
  ! CH4 : 0.85 * qsat because of dissolution in N2
  qsatCH4 = 0.85 * qsatCH4

  ! Flux at the surface [kg.m-2.s-1] :
  flux = rhoair * Cd * ws

  ! <North Pole>
  if (REAL(latitude_deg(ig)) .ge. 70.) then
    flux = Flnp * flux
  ! <South Pole>
  else if (REAL(latitude_deg(ig)) .le. -70.) then
    flux = Flsp * flux
  ! <Mid Latitudes>
  else
    flux = 0.0
  endif

  ! Empty tank ?
  if (tankCH4(ig) .le. 1.e-30) then
    flux = 0.0
    tankCH4(ig) = 1.e-30
  endif
  
  ! Flux of CH4 at the surface [kg.m-2.s-1.mol.mol-1] :
  fluxCH4 = flux * (qsatCH4 - pqCH4(ig,1))


!------------------------------------
! 2. IMPLICIT SCHEME
!------------------------------------

  ! 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. 0.) then
    tankCH4(ig) = tankCH4(ig) + dtankCH4
  else
    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. - tsurf(ig) / 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 * 190.4 * (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

  ! >>> [TEMPO : BBT]
  !open(501,file='Evap_CH4.out')
  !if(REAL(latitude_deg(ig)) .ge. 70.) then
  !      write(501,*) '-----------------------------------'
  !      write(501,*) 'Latitude =', REAL(latitude_deg(ig))
  !      write(501,*) ''
  !      write(501,*) 'Initialisation variables :'
  !     write(501,*) 'Air density =', rhoair, 'kg.m-3'
  !      write(501,*) 'Horizontal wind =', ws, 'm.s-1'
  !      write(501,*) 'Turbulent term =', Cd
  !      write(501,*) 'CH4 saturation at the surface =', qsatCH4, 'mol.mol-1'
  !      write(501,*) ''
  !      write(501,*) 'Flux :'
  !      write(501,*) 'Flux of CH4 =', fluxCH4, 'kg.m-2.s-1.mol.mol-1'
  !      write(501,*) ''
  !      write(501,*) 'Variation of CH4 :'
  !      write(501,*) 'New molar fraction of CH4 =', newpqCH4, 'mol.mol-1'
  !      write(501,*) 'Trend of CH4 =', pdqCH4(ig), 'mol.mol-1.s-1'
  !      write(501,*) 'Variation of CH4 tank =', dtankCH4, 'm'
  !      write(501,*) ''
  !      write(501,*) 'Latent Heat :'
  !      write(501,*) 'Latente heat of CH4 vaporisation =', LvCH4, 'J.kg-1.mol.mol-1.s-1'
  !      write(501,*) 'dtsurfevap :', dtsurfevap(ig), 'K.s-1'
  !      write(501,*) '-----------------------------------'
  !endif
  !close(501)
  ! <<< [TEMPO : BBT]

enddo ! End of main loop on the horizontal grid

return

end subroutine evapCH4