source: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 @ 5551

MODULE lmdz_aviation
!----------------------------------------------------------------
! Module for aviation and contrails

IMPLICIT NONE

CONTAINS

SUBROUTINE aviation_water_emissions( &
      klon, klev, dtime, pplay, temp, qtot, &
      flight_h2o, d_q_avi &
      )

USE lmdz_lscp_ini, ONLY: RD

IMPLICIT NONE

INTEGER,                      INTENT(IN)  :: klon, klev ! number of horizontal grid points and vertical levels
REAL,                         INTENT(IN)  :: dtime      ! time step [s]
REAL, DIMENSION(klon,klev),   INTENT(IN)  :: pplay      ! mid-layer pressure [Pa]
REAL, DIMENSION(klon,klev),   INTENT(IN)  :: temp       ! temperature (K)
REAL, DIMENSION(klon,klev),   INTENT(IN)  :: qtot       ! total specific humidity (in vapor phase) [kg/kg] 
REAL, DIMENSION(klon,klev),   INTENT(IN)  :: flight_h2o ! aviation emitted H2O concentration [kgH2O/s/m3]
REAL, DIMENSION(klon,klev),   INTENT(OUT) :: d_q_avi    ! water vapor tendency from aviation [kg/kg]
! Local
INTEGER :: i, k
REAL :: rho

DO i=1, klon
  DO k=1, klev
    !--Dry density [kg/m3]
    rho = pplay(i,k) / temp(i,k) / RD
    !--Dry air mass [kg/m2]
    !rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
    !--Cell thickness [m]
    !dz = rhodz / rho
    !--Cell dry air mass [kg]
    !M_cell = rhodz * cell_area(i)

    !--q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
    ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )  
    !      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) 
    !--The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) 
    !--flight_h2O is in kg H2O / s / m3
    
    !d_q_avi(i,k) = ( M_cell * qtot(i,k) + V_cell * flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
    !             / ( M_cell             + V_cell * flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
    !             - qtot(i,k)
    !--NB., M_cell = V_cell * rho
    !d_q_avi(i,k) = ( rho * qtot(i,k) + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
    !             / ( rho             + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
    !             - qtot(i,k)
    !--Same formula, more computationally effective but less readable
    d_q_avi(i,k) = flight_h2o(i,k) * ( 1. - qtot(i,k) ) &
                 / ( rho / dtime / ( 1. - qtot(i,k) ) + flight_h2o(i,k) )
  ENDDO
ENDDO

END SUBROUTINE aviation_water_emissions


!**********************************************************************************
SUBROUTINE contrails_formation_evolution( &
      dtime, pplay, temp, qsat, qsatl, gamma_cond, rcont_seri, flight_dist, &
      cldfra, qvc, dz, pdf_loc, pdf_scale, pdf_alpha, &
      Tcritcont, qcritcont, potcontfraP, potcontfraNP, contfra, &
      dcontfra_cir, dcf_avi, dqvc_avi, dqi_avi &
      )

USE lmdz_lscp_ini, ONLY: RCPD, EPS_W, RTT
USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation
USE lmdz_lscp_ini, ONLY: linear_contrails_lifetime
USE lmdz_lscp_ini, ONLY: eps

USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf

IMPLICIT NONE

!
! Input
!
REAL, INTENT(IN)  :: dtime        ! time step [s]
REAL, INTENT(IN)  :: pplay        ! layer pressure [Pa]
REAL, INTENT(IN)  :: temp         ! temperature [K]
REAL, INTENT(IN)  :: qsat         ! saturation specific humidity [kg/kg]
REAL, INTENT(IN)  :: qsatl        ! saturation specific humidity w.r.t. liquid [kg/kg]
REAL, INTENT(IN)  :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
REAL, INTENT(IN)  :: rcont_seri   ! ratio of contrails fraction to total cloud fraction [-]
REAL, INTENT(IN)  :: flight_dist  ! aviation distance flown concentration [m/s/m3]
REAL, INTENT(IN)  :: cldfra       ! cloud fraction [-]
REAL, INTENT(IN)  :: qvc          ! gridbox-mean vapor in the cloud [kg/kg]
REAL, INTENT(IN)  :: dz           ! cell width [m]
REAL, INTENT(IN)  :: pdf_loc      ! location parameter of the clear sky PDF [%]
REAL, INTENT(IN)  :: pdf_scale    ! scale parameter of the clear sky PDF [%]
REAL, INTENT(IN)  :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
!
! Output
!
REAL, INTENT(OUT) :: Tcritcont    ! critical temperature for contrail formation [K]
REAL, INTENT(OUT) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
REAL, INTENT(OUT) :: potcontfraP  ! potential persistent contrail fraction [-]
REAL, INTENT(OUT) :: potcontfraNP ! potential non-persistent contrail fraction [-]
REAL, INTENT(OUT) :: contfra      ! linear contrail fraction [-]
REAL, INTENT(OUT) :: dcontfra_cir ! linear contrail fraction to cirrus cloud fraction tendency [s-1]
REAL, INTENT(OUT) :: dcf_avi      ! cloud fraction tendency because of aviation [s-1]
REAL, INTENT(OUT) :: dqvc_avi     ! specific ice content tendency because of aviation [kg/kg/s]
REAL, INTENT(OUT) :: dqi_avi      ! specific cloud water vapor tendency because of aviation [kg/kg/s]
!
! Local
!
! for Schmidt-Appleman criteria
REAL, DIMENSION(1) :: ztemp, zpplay, qzero, zqsatl, zdqsatl
REAL :: Gcont, qsatl_crit, psatl_crit, pcrit
REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3
REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc
REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat, pdf_q_above_qnuc
REAL :: qpotcontP
!
! for new contrail formation
REAL :: contrail_cross_section, contfra_new

qzero(:) = 0.

!---------------------------------
!--  SCHIMDT-APPLEMAN CRITERIA  --
!---------------------------------
!--Revised by Schumann (1996) and Rap et al. (2010)

!--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute 
!--in Pa / K. See Rap et al. (2010) in JGR.
!--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
Gcont = EI_H2O_aviation * RCPD * pplay &
       / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
!--critical temperature below which no liquid contrail can form in exhaust
!--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2
!--in Kelvins
Tcritcont = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
       + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2

ztemp(1) = Tcritcont
zpplay(1) = pplay
CALL calc_qsat_ecmwf(1, ztemp, qzero, zpplay, RTT, 1, .FALSE., zqsatl, zdqsatl)
qsatl_crit = zqsatl(1)

psatl_crit = qsatl_crit / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit ) * pplay
pcrit = Gcont * ( temp - Tcritcont ) + psatl_crit
qcritcont = EPS_W * pcrit / ( pplay - ( 1. - EPS_W ) * pcrit )


IF ( ( ( 1. - cldfra ) .GT. eps ) .AND. ( temp .LT. Tcritcont ) ) THEN

  pdf_x = qcritcont / qsatl * 100.
  pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
  pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
  pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
  pdf_fra_above_qcritcont = EXP( - pdf_y ) * ( 1. - cldfra )
  pdf_q_above_qcritcont = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qcritcont ) * qsatl / 100.

  pdf_x = gamma_cond * qsat / qsatl * 100.
  pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
  pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
  pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
  pdf_fra_above_qnuc = EXP( - pdf_y ) * ( 1. - cldfra )
  pdf_q_above_qnuc = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qnuc ) * qsatl / 100.

  pdf_x = qsat / qsatl * 100.
  pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
  pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
  pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
  pdf_fra_above_qsat = EXP( - pdf_y ) * ( 1. - cldfra )
  pdf_q_above_qsat = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qsat ) * qsatl / 100.

  potcontfraNP = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
  potcontfraP = MAX(0., MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
                    pdf_fra_above_qcritcont - pdf_fra_above_qnuc))
  qpotcontP = MAX(0., MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
                  pdf_q_above_qcritcont - pdf_q_above_qnuc))

ELSE

  potcontfraNP = 0.
  potcontfraP = 0.

ENDIF ! temp .LT. Tcritcont

!--Convert existing contrail fraction into "natural" cirrus cloud fraction
contfra = rcont_seri * cldfra
dcontfra_cir = contfra * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
contfra = contfra + dcontfra_cir

!--Add a source of contrails from aviation
dcf_avi = 0.
dqi_avi = 0.
dqvc_avi = 0.
IF ( potcontfraP .GT. eps ) THEN
  contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA(dz)
  contfra_new = MIN(1., flight_dist * dtime * contrail_cross_section)
  dcf_avi = potcontfraP * contfra_new
  IF ( cldfra .GT. eps ) THEN
    dqvc_avi = dcf_avi * qvc / cldfra
  ELSE
    dqvc_avi = dcf_avi * qsat
  ENDIF
  dqi_avi = dcf_avi * qpotcontP / potcontfraP - dqvc_avi

  !--Add tendencies
  contfra = contfra + dcf_avi
ENDIF

END SUBROUTINE contrails_formation_evolution


!**********************************************************************************
SUBROUTINE contrails_mixing( &
      dtime, pplay, shear, pbl_eps, cell_area, dz, temp, qtot, qsat, &
      subfra, qsub, issrfra, qissr, cldfra, qcld, qvc, rcont_seri, &
      dcf_mix, dqvc_mix, dqi_mix, dqt_mix, dcf_mix_issr, dqt_mix_issr &
      )

!----------------------------------------------------------------------
! This subroutine calculates the mixing of contrails in their environment.
! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
! December 2024
!----------------------------------------------------------------------

USE lmdz_lscp_ini,        ONLY: RPI, eps, ok_unadjusted_clouds
USE lmdz_lscp_ini,        ONLY: aspect_ratio_contrails, coef_mixing_contrails, &
                                coef_shear_contrails, chi_mixing_contrails

IMPLICIT NONE

!
! Input
!
REAL, INTENT(IN)    :: dtime        ! time step [s]
!
REAL, INTENT(IN)    :: pplay        ! layer pressure [Pa]
REAL, INTENT(IN)    :: shear        ! vertical shear [s-1]
REAL, INTENT(IN)    :: pbl_eps      ! TKE dissipation[m2/s3]
REAL, INTENT(IN)    :: cell_area    ! cell area [m2]
REAL, INTENT(IN)    :: dz           ! cell width [m]
REAL, INTENT(IN)    :: temp         ! temperature [K]
REAL, INTENT(IN)    :: qtot         ! total specific humidity (without precip) [kg/kg]
REAL, INTENT(IN)    :: qsat         ! saturation specific humidity [kg/kg]
REAL, INTENT(IN)    :: subfra       ! subsaturated clear sky fraction [-]
REAL, INTENT(IN)    :: qsub         ! gridbox-mean subsaturated clear sky specific water [kg/kg]
REAL, INTENT(IN)    :: issrfra      ! ISSR fraction [-]
REAL, INTENT(IN)    :: qissr        ! gridbox-mean ISSR specific water [kg/kg]
REAL, INTENT(IN)    :: cldfra       ! cloud fraction [-]
REAL, INTENT(IN)    :: qcld         ! gridbox-mean cloudy specific total water [kg/kg]
REAL, INTENT(IN)    :: qvc          ! gridbox-mean cloudy specific water vapor [kg/kg]
REAL, INTENT(IN)    :: rcont_seri   ! ratio of contrails fraction to total cloud fraction [-]
!
!  Input/Output
!
REAL, INTENT(INOUT) :: dcf_mix      ! cloud fraction tendency because of cloud mixing [s-1]
REAL, INTENT(INOUT) :: dqvc_mix     ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
REAL, INTENT(INOUT) :: dqi_mix      ! specific ice content tendency because of cloud mixing [kg/kg/s]
REAL, INTENT(INOUT) :: dqt_mix      ! total water tendency because of cloud mixing [kg/kg/s]
REAL, INTENT(INOUT) :: dcf_mix_issr ! cloud fraction tendency because of cloud mixing in ISSR [s-1]
REAL, INTENT(INOUT) :: dqt_mix_issr ! total water tendency because of cloud mixing in ISSR [kg/kg/s]
!
! Local
!
REAL :: dqt_mix_sub_cont, dqt_mix_issr_cont
REAL :: dcf_mix_sub_cont, dcf_mix_issr_cont
REAL :: dqvc_mix_sub_cont, dqvc_mix_issr_cont
REAL :: dcf_mix_cont, dqvc_mix_cont, dqi_mix_cont, dqt_mix_cont
REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
REAL :: envfra_mix, cldfra_mix
REAL :: L_shear, shear_fra
REAL :: sigma_mix, issrfra_mix, subfra_mix
REAL :: qvapincld, qvapinmix, qvapincld_new, qiceincld

!--Initialisation
dcf_mix_sub_cont   = 0.
dqt_mix_sub_cont   = 0.
dqvc_mix_sub_cont  = 0.
dcf_mix_issr_cont  = 0.
dqt_mix_issr_cont  = 0.
dqvc_mix_issr_cont = 0.

!-- PART 1 - TURBULENT DIFFUSION

!--Clouds within the mesh are assumed to be ellipses. The length of the
!--semi-major axis is a and the length of the semi-minor axis is b.
!--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
!--In particular, it is 0 if cldfra = 1.
!--clouds_perim is the total perimeter of the clouds within the mesh,
!--not considering interfaces with other meshes (only the interfaces with clear
!--sky are taken into account).
!--
!--The area of each cloud is A = a * b * RPI,
!--and the perimeter of each cloud is
!-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
!--
!--With cell_area the area of the cell, we have:
!-- cldfra = A * N_cld_mix / cell_area
!-- clouds_perim = P * N_cld_mix
!--
!--We assume that the ratio between b and a is a function of
!--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
!--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
!--are spherical.
!-- b / a = bovera = MAX(0.1, cldfra)
bovera = aspect_ratio_contrails
!--P / a is a function of b / a only, that we can calculate
!-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
!--The clouds perimeter is imposed using the formula from Morcrette 2012,
!--based on observations.
!-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
!--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
!-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
!-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - cldfra ) / bovera
!--and finally,
!-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
N_cld_mix = coef_mixing_contrails * cldfra * ( 1. - cldfra ) * cell_area &
          / Povera / a_mix

!--The time required for turbulent diffusion to homogenize a region of size
!--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
!--We compute L_mix and assume that the cloud is mixed over this length
L_mix = SQRT( dtime**3 * pbl_eps )
!--The mixing length cannot be greater than the semi-minor axis. In this case,
!--the entire cloud is mixed.
L_mix = MIN(L_mix, a_mix * bovera)

!--The fraction of clear sky mixed is
!-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
envfra_mix = N_cld_mix * RPI / cell_area &
           * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
!--The fraction of cloudy sky mixed is
!-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
cldfra_mix = N_cld_mix * RPI / cell_area &
           * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )


!-- PART 2 - SHEARING

!--The clouds are then sheared. We keep the shape and number
!--assumptions from before. The clouds are sheared with a random orientation
!--of the wind, on average we assume that the wind and the semi-major axis
!--make a 45 degrees angle. Moreover, the contrails only mix
!--along their semi-minor axis (b), because it is easier to compute.
!--With this, the clouds increase in size along b only, by a factor
!--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
L_shear = coef_shear_contrails * shear * dz * dtime
!--therefore, the fraction of clear sky mixed is
!-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
!--and the fraction of cloud mixed is
!-- N_cld_mix * ( a * b - a * (b - L_shear * SQRT(2.) / 2.) ) * RPI / 2. / cell_area
!--(note that they are equal)
shear_fra = RPI * SQRT(2.) / 2. * L_shear * a_mix / 2. * N_cld_mix / cell_area
!--and the environment and cloud mixed fractions are the same,
!--which we add to the previous calculated mixed fractions.
!--We therefore assume that the sheared clouds and the turbulent
!--mixed clouds are different.
envfra_mix = envfra_mix + shear_fra
cldfra_mix = cldfra_mix + shear_fra


!-- PART 3 - CALCULATION OF THE MIXING PROPERTIES

!--The environment fraction is allocated to subsaturated sky or supersaturated sky,
!--according to the factor sigma_mix. This is computed as the ratio of the
!--subsaturated sky fraction to the environment fraction, corrected by a factor
!--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
!--supersaturated sky is favoured. Physically, this means that it is more likely
!--to have supersaturated sky around the cloud than subsaturated sky.
sigma_mix = subfra / ( subfra + chi_mixing_contrails * issrfra )
subfra_mix = MIN( sigma_mix * envfra_mix, subfra )
issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra )
cldfra_mix = MIN( cldfra_mix, cldfra )

!--First, we mix the subsaturated sky (subfra_mix) and the cloud close
!--to this fraction (sigma_mix * cldfra_mix).
IF ( subfra .GT. eps ) THEN
  !--We compute the total humidity in the mixed air, which
  !--can be either sub- or supersaturated.
  qvapinmix = ( qsub * subfra_mix / subfra &
            + qcld * cldfra_mix * sigma_mix / cldfra ) &
            / ( subfra_mix + cldfra_mix * sigma_mix )

  IF ( ok_unadjusted_clouds ) THEN
    qiceincld = ( qcld - qvc ) * cldfra_mix * sigma_mix / cldfra &
              / ( subfra_mix + cldfra_mix * sigma_mix )
    qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
        qvapinmix, qiceincld, temp, qsat, pplay, dtime)
    IF ( qvapincld_new .EQ. 0. ) THEN
      !--If all the ice has been sublimated, we sublimate
      !--completely the cloud and do not activate the sublimation
      !--process
      !--Tendencies and diagnostics
      dcf_mix_sub_cont = - sigma_mix * cldfra_mix
      dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
      dqvc_mix_sub_cont = dcf_mix_sub_cont * qvc / cldfra
    ELSE
      dcf_mix_sub_cont = subfra_mix
      dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
      dqvc_mix_sub_cont = dcf_mix_sub_cont * qvapincld_new
    ENDIF
  ELSE
    IF ( qvapinmix .GT. qsat ) THEN
      !--If the mixed air is supersaturated, we condense the subsaturated
      !--region which was mixed.
      dcf_mix_sub_cont = subfra_mix
      dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
      dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
    ELSE
      !--Else, we sublimate the cloud which was mixed.
      dcf_mix_sub_cont = - sigma_mix * cldfra_mix
      dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
      dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
    ENDIF
  ENDIF ! ok_unadjusted_clouds
ENDIF ! subfra .GT. eps

!--We then mix the supersaturated sky (issrfra_mix) and the cloud,
!--for which the mixed air is always supersatured, therefore
!--the cloud necessarily expands
IF ( issrfra .GT. eps ) THEN

  IF ( ok_unadjusted_clouds ) THEN
    qvapinmix = ( qissr * issrfra_mix / issrfra &
              + qcld * cldfra_mix * (1. - sigma_mix) / cldfra ) &
              / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
    qiceincld = ( qcld - qvc ) * cldfra_mix * (1. - sigma_mix) / cldfra &
              / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
    qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
        qvapinmix, qiceincld, temp, qsat, pplay, dtime)
    dcf_mix_issr_cont = issrfra_mix
    dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
    dqvc_mix_issr_cont = dcf_mix_issr_cont * qvapincld_new
  ELSE
    !--In this case, the additionnal vapor condenses
    dcf_mix_issr_cont = issrfra_mix
    dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
    dqvc_mix_issr_cont = dcf_mix_issr_cont * qsat
  ENDIF ! ok_unadjusted_clouds

ENDIF ! issrfra .GT. eps

!--Sum up the tendencies from subsaturated sky and supersaturated sky
dcf_mix_cont  = dcf_mix_sub_cont  + dcf_mix_issr_cont
dqt_mix_cont  = dqt_mix_sub_cont  + dqt_mix_issr_cont
dqvc_mix_cont = dqvc_mix_sub_cont + dqvc_mix_issr_cont
dqi_mix_cont  = dqt_mix_cont - dqvc_mix_cont

!--Outputs
!--The mixing tendencies are a linear combination of the calculation done for "classical" cirrus
!--and contrails
dcf_mix = ( 1. - rcont_seri ) * dcf_mix + rcont_seri * dcf_mix_cont
dqvc_mix = ( 1. - rcont_seri ) * dqvc_mix + rcont_seri * dqvc_mix_cont
dqi_mix = ( 1. - rcont_seri ) * dqi_mix + rcont_seri * dqi_mix_cont
dqt_mix = ( 1. - rcont_seri ) * dqt_mix + rcont_seri * dqt_mix_cont
dcf_mix_issr = ( 1. - rcont_seri ) * dcf_mix_issr + rcont_seri * dcf_mix_issr_cont
dqt_mix_issr = ( 1. - rcont_seri ) * dqt_mix_issr + rcont_seri * dqt_mix_issr_cont

END SUBROUTINE contrails_mixing


!**********************************************************************************
FUNCTION qvapincld_depsub_contrails( &
    qvapincld, qiceincld, temp, qsat, pplay, dtime)

USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, EPS_W
USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
USE lmdz_lscp_ini,        ONLY: rm_ice_crystals_contrails

IMPLICIT NONE

! inputs
REAL :: qvapincld     ! 
REAL :: qiceincld     ! 
REAL :: temp          ! 
REAL :: qsat          ! 
REAL :: pplay         ! 
REAL :: dtime         ! time step [s]
! output
REAL :: qvapincld_depsub_contrails
! local
REAL :: qice_ratio
REAL :: pres_sat, rho, kappa
REAL :: air_thermal_conduct, water_vapor_diff
REAL :: rm_ice

!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
!
!--The deposition equation is
!-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
!--from Lohmann et al. (2016), where
!--alpha is the deposition coefficient [-]
!--mi is the mass of one ice crystal [kg]
!--C is the capacitance of an ice crystal [m]
!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
!--R_v is the specific gas constant for humid air [J/kg/K]
!--T is the temperature [K]
!--esi is the saturation pressure w.r.t. ice [Pa]
!--Dv is the diffusivity of water vapor [m2/s]
!--Ls is the specific latent heat of sublimation [J/kg/K]
!--ka is the thermal conductivity of dry air [J/m/s/K]
!
!--alpha is a coefficient to take into account the fact that during deposition, a water
!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
!--coherent (with the same structure). It has no impact for sublimation.
!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
!--during deposition, and alpha = 1. during sublimation.
!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
!-- C = capa_cond_cirrus * rm_ice
!
!--We have qice = Nice * mi, where Nice is the ice crystal
!--number concentration per kg of moist air
!--HYPOTHESIS 1: the ice crystals are spherical, therefore
!-- mi = 4/3 * pi * rm_ice**3 * rho_ice
!--HYPOTHESIS 2: the ice crystals are monodisperse with the
!--initial radius rm_ice_0.
!--NB. this is notably different than the assumption
!--of a distributed qice in the cloud made in the sublimation process;
!--should it be consistent?
!
!--As the deposition process does not create new ice crystals,
!--and because we assume a same rm_ice value for all crystals
!--therefore the sublimation process does not destroy ice crystals
!--(or, in a limit case, it destroys all ice crystals), then
!--Nice is a constant during the sublimation/deposition process.
!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
!
!--The deposition equation then reads:
!-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
!-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
!--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
!--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
!-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
!--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
!--and we have
!-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
!-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
!--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
!
!--This system of equations can be resolved with an exact
!--explicit numerical integration, having one variable resolved
!--explicitly, the other exactly. The exactly resolved variable is
!--the one decreasing, so it is qvc if the process is
!--condensation, qi if it is sublimation.
!
!--kappa is computed as an initialisation constant, as it depends only
!--on temperature and other pre-computed values
pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
      * ( RV * temp / water_vapor_diff / pres_sat &
        + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process

!--Here, the initial vapor in the cloud is qvapincld, and we compute
!--the new vapor qvapincld_depsub_contrails

rm_ice = rm_ice_crystals_contrails

IF ( qvapincld .GE. qsat ) THEN
  !--If the cloud is initially supersaturated
  !--Exact explicit integration (qvc exact, qice explicit)
  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
                * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
ELSE
  !--If the cloud is initially subsaturated
  !--Exact explicit integration (qice exact, qvc explicit)
  !--The barrier is set so that the resulting vapor in cloud
  !--cannot be greater than qsat
  !--qice_ratio is the ratio between the new ice content and
  !--the old one, it is comprised between 0 and 1
  qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )

  IF ( qice_ratio .LT. 0. ) THEN
    !--The new vapor in cloud is set to 0 so that the
    !--sublimation process does not activate
    qvapincld_depsub_contrails = 0.
  ELSE
    !--Else, the sublimation process is activated with the
    !--diagnosed new cloud water vapor
    !--The new vapor in the cloud is increased with the
    !--sublimated ice
    qvapincld_depsub_contrails = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
    !--The new vapor in the cloud cannot be greater than qsat
    qvapincld_depsub_contrails = MIN(qvapincld_depsub_contrails, qsat)
  ENDIF ! qice_ratio .LT. 0.
ENDIF ! qvapincld .GT. qsat

END FUNCTION qvapincld_depsub_contrails


!**********************************************************************************
FUNCTION contrail_cross_section_onera(dz)

USE lmdz_lscp_ini, ONLY: initial_width_contrails

IMPLICIT NONE

!
! Input
!
REAL :: dz ! cell width [m]
!
! Output
!
REAL :: contrail_cross_section_onera ! [m2]
!
! Local
!

contrail_cross_section_onera = initial_width_contrails * dz

END FUNCTION contrail_cross_section_onera

SUBROUTINE read_aviation_emissions(klon, klev, flight_dist, flight_h2o)
    ! This subroutine allows to read the traffic density data read in the file aviation.nc
    ! This file is defined in ./COMP/lmdz.card
    ! Field names in context_input_lmdz.xml should be the same as in the file.
    USE netcdf
    USE mod_phys_lmdz_para
    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, nbp_lev_=>nbp_lev, &
                                 klon_glo, grid2Dto1D_glo, grid_type, unstructured
    USE lmdz_xios
    USE print_control_mod, ONLY: lunout
    IMPLICIT NONE

    INTEGER,                    INTENT(IN)  :: klon, klev  ! number of horizontal grid points and vertical levels
    REAL, DIMENSION(klon,klev), INTENT(OUT) :: flight_dist ! Aviation distance flown concentration [m/s/m3]
    REAL, DIMENSION(klon,klev), INTENT(OUT) :: flight_h2o  ! Aviation emitted H2O [kgH2O/s/m3]

    !----------------------------------------------------
    ! Local variable
    !----------------------------------------------------
    REAL, DIMENSION(klon_mpi,klev,1) :: flight_dist_mpi

    !--Initialisation
    flight_dist(:,:) = 0.
    flight_h2o(:,:) = 0.

    ! Read the data from the file
    ! is_omp_master is necessary to make XIOS works
    IF (is_omp_master) CALL xios_recv_field("KMFLOWN_interp", flight_dist_mpi(:,:,1))

    ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev)
    ! (klon_mpi,klon) = (200,50) avec 80 MPI, 4 OMP, nbp40
    CALL scatter_omp(flight_dist_mpi(:,:,1), flight_dist)

END SUBROUTINE read_aviation_emissions

END MODULE lmdz_aviation