Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5453)
@@ -7,9 +7,59 @@
 CONTAINS
 
+SUBROUTINE aviation_water_emissions( &
+      klon, klev, dtime, paprs, pplay, temp, qtot, cell_area, &
+      flight_h2o, d_q_avi &
+      )
+
+USE lmdz_lscp_ini, ONLY: RD, RG
+
+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+1), INTENT(IN)  :: paprs      ! inter-layer pressure [Pa]
+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),        INTENT(IN)  :: cell_area  ! area of each cell [m2]
+REAL, DIMENSION(klon,klev),   INTENT(IN)  :: flight_h2o ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
+REAL, DIMENSION(klon,klev),   INTENT(OUT) :: d_q_avi    ! water vapor tendency from aviation [kg/kg]
+! Local
+INTEGER :: i, k
+REAL :: rho, rhodz, dz, M_cell
+
+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 / mesh
+    
+    !d_q_avi(i,k) = ( M_cell * qtot(i,k) + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
+    !             / ( M_cell             + 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) ) &
+                 / ( M_cell / 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, V_cell, M_cell, pdf_loc, pdf_scale, pdf_alpha, &
+      cldfra, qvc, dz, V_cell, pdf_loc, pdf_scale, pdf_alpha, &
       Tcritcont, qcritcont, potcontfraP, potcontfraNP, contfra, &
       dcf_avi, dqvc_avi, dqi_avi &
@@ -35,9 +85,9 @@
 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  ! 
+REAL, INTENT(IN)  :: flight_dist  ! aviation distance flown within the mesh [m/s/mesh]
 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)  :: V_cell       ! cell volume [m3]
-REAL, INTENT(IN)  :: M_cell       ! cell mass [kg]
 REAL, INTENT(IN)  :: pdf_loc      ! location parameter of the clear sky PDF [%]
 REAL, INTENT(IN)  :: pdf_scale    ! scale parameter of the clear sky PDF [%]
@@ -46,12 +96,12 @@
 ! Output
 !
-REAL, INTENT(OUT) :: Tcritcont    !
-REAL, INTENT(OUT) :: qcritcont    !
-REAL, INTENT(OUT) :: potcontfraP  !
-REAL, INTENT(OUT) :: potcontfraNP !
-REAL, INTENT(OUT) :: contfra      !
-REAL, INTENT(OUT) :: dcf_avi      !
-REAL, INTENT(OUT) :: dqvc_avi     !
-REAL, INTENT(OUT) :: dqi_avi      !
+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      ! contrail fraction [-]
+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
@@ -65,19 +115,8 @@
 REAL :: qpotcontP
 !
-!
+! for new contrail formation
 REAL :: contrail_cross_section, contfra_new
 
 qzero(:) = 0.
-
-!--more local variables for diagnostics
-!--values from Schumann, Meteorol Zeitschrift, 1996
-!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
-!--Qheat = 43.  /  50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
-!REAL, PARAMETER :: EiH2O=1.25  !--emission index of water vapour for kerosene (kg kg-1)
-!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) 
-!REAL, PARAMETER :: eta=0.3     !--average propulsion efficiency of the aircraft
-!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute 
-!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
-!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) 
 
 !---------------------------------
@@ -86,9 +125,11 @@
 !--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)
-!--in Pa / K
 Gcont = EI_H2O_aviation * RCPD * pplay &
        / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
-!--critical T_LM below which no liquid contrail can form in exhaust
+!--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 ) &
@@ -105,5 +146,5 @@
 
 
-IF ( temp .LT. Tcritcont ) THEN !--contrail formation is possible
+IF ( ( ( 1. - cldfra ) .GT. eps ) .AND. ( temp .LT. Tcritcont ) ) THEN
 
   pdf_x = qcritcont / qsatl * 100.
@@ -132,8 +173,8 @@
 
   potcontfraNP = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
-  potcontfraP = MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
-                    pdf_fra_above_qcritcont - pdf_fra_above_qnuc)
-  qpotcontP = MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
-                  pdf_q_above_qcritcont - pdf_q_above_qnuc)
+  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
@@ -143,16 +184,4 @@
 
 ENDIF ! temp .LT. Tcritcont
-
-
-! Ajout des émissions de H2O dues à l'aviation
-! 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 / cell
-!
-!IF (ok_plane_h2o) THEN 
-!  q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) 
-!ENDIF
 
 !--Convert existing contrail fraction into "natural" cirrus cloud fraction
@@ -165,5 +194,5 @@
 dqvc_avi = 0.
 IF ( potcontfraP .GT. eps ) THEN
-  contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
+  contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA(dz)
   contfra_new = MIN(1., flight_dist * dtime * contrail_cross_section / V_cell)
   dcf_avi = potcontfraP * contfra_new
@@ -565,16 +594,62 @@
 
 !**********************************************************************************
-FUNCTION contrail_cross_section_onera()
+FUNCTION contrail_cross_section_onera(dz)
+
+USE lmdz_lscp_init, ONLY: initial_width_contrails
 
 IMPLICIT NONE
 
-! input
-! output
+!
+! Input
+!
+REAL :: dz ! cell width [m]
+!
+! Output
+!
 REAL :: contrail_cross_section_onera ! [m2]
-! local
-
-contrail_cross_section_onera = 200. * 200.
+!
+! Local
+!
+
+contrail_cross_section_onera = initial_width_contrails * dz
 
 END FUNCTION contrail_cross_section_onera
+
+
+SUBROUTINE read_aviation_emissions( &
+        klon, klev, latitude_deg, longitude_deg, pplay, &
+        flight_dist, flight_h2o &
+        )
+
+IMPLICIT NONE
+!
+! Input
+!
+INTEGER, INTENT(IN)                        :: klon, klev       ! number of horizontal grid points and vertical levels
+REAL,    INTENT(IN),  DIMENSION(klon)      :: latitude_deg     ! latitude of the grid points [deg]
+REAL,    INTENT(IN),  DIMENSION(klon)      :: longitude_deg    ! longitude of the grid points [deg]
+REAL,    INTENT(IN),  DIMENSION(klon,klev) :: pplay            ! layer pressure [Pa]
+!
+! Output
+!
+REAL,    INTENT(OUT), DIMENSION(klon,klev) :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
+REAL,    INTENT(OUT), DIMENSION(klon,klev) :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
+!
+! Local
+!
+INTEGER :: i, k
+
+!--Initialisation
+flight_dist(:,:) = 0.
+flight_h2o(:,:) = 0.
+
+!DO i=1, klon
+! IF ( ( latitude_deg(i) .GE. 42. ) .AND. ( latitude_deg(i) .LE. 48. ) ) THEN 
+!   flight_dist(i,14) = 50000.  !--5000 m of flight/second in grid cell x 10 scaling
+!   flight_h2o(i,14) = 100.  !--10 kgH2O/second in grid cell x 10 scaling
+! ENDIF
+!ENDDO
+
+END SUBROUTINE read_aviation_emissions
 
 END MODULE lmdz_aviation
@@ -754,43 +829,2 @@
 !!
 !END SUBROUTINE airplane
-!
-!!********************************************************************
-!! simple routine to initialise flight_m and test a flight corridor
-!!--Olivier Boucher - 2021 
-!!
-!SUBROUTINE flight_init()
-!  USE dimphy
-!  USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
-!  IMPLICIT NONE
-!  INTEGER :: i
-!
-!  ALLOCATE(flight_m(klon,klev))
-!  ALLOCATE(flight_h2o(klon,klev))
-!  !
-!  flight_m(:,:) = 0.0    !--initialisation
-!  flight_h2o(:,:) = 0.0  !--initialisation
-!  !
-!  DO i=1, klon
-!   IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN 
-!     flight_m(i,38) = 50000.0  !--5000 m of flight/second in grid cell x 10 scaling
-!   ENDIF
-!  ENDDO
-!  
-!  RETURN
-!END SUBROUTINE flight_init
-!
-!     !--add a source of cirrus from aviation contrails
-!     IF (ok_plane_contrail) THEN
-!       drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell    !--tendency rneb due to aviation [s-1]
-!       drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime)                     !--majoration
-!       dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss)                          !--tendency q aviation [kg kg-1 s-1]
-!       rneb = rneb + drneb_avi(i,k)*dtime                                     !--add tendency to rneb
-!       qcld = qcld + dqss_avi*dtime                                           !--add tendency to qcld
-!       rnebss = rnebss - drneb_avi(i,k)*dtime                                 !--add tendency to rnebss
-!       qss = qss - dqss_avi*dtime                                             !--add tendency to qss
-!     ELSE 
-!       drneb_avi(i,k)=0.0 
-!     ENDIF
-!
-!  RETURN
-!END SUBROUTINE ice_sursat
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5453)
@@ -175,6 +175,6 @@
   !--------------------------------------------------
   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rcont_seri       ! ratio of contrails fraction to total cloud fraction [-]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
   
   ! OUTPUT variables
@@ -1066,5 +1066,5 @@
           rvc_seri(i,k) = qvc(i) / zq(i)
         ELSE
-          rvc_seri(i,k) = min_ratio
+          rvc_seri(i,k) = 0.
         ENDIF
         !--The MIN barrier is NEEDED because of:
@@ -1077,5 +1077,5 @@
             rcont_seri(i,k) = contfra(i) / rneb(i,k)
           ELSE
-            rcont_seri(i,k) = min_ratio
+            rcont_seri(i,k) = 0.
           ENDIF
           !--This barrier should never be activated
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5453)
@@ -157,6 +157,6 @@
 !
 REAL,     INTENT(INOUT), DIMENSION(klon) :: rcont_seri    ! ratio of contrails fraction to total cloud fraction [-]
-REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   ! 
-REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    ! 
+REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   ! aviation distance flown within the mesh [m/s/mesh]
+REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
 !
 !  Output
@@ -238,6 +238,5 @@
 !
 ! for cell properties
-REAL :: rho, rhodz, dz
-REAL :: V_cell, M_cell
+REAL :: rho, rhodz, dz, V_cell
 
 qzero(:) = 0.
@@ -336,6 +335,4 @@
       !--Cell volume [m3]
       V_cell = dz * cell_area(i)
-      !--Cell dry air mass [kg]
-      M_cell = rhodz * cell_area(i)
 
 
@@ -504,6 +501,5 @@
         ENDIF
         pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
-        pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
-        pdf_alpha = MIN(10., pdf_alpha)
+        pdf_alpha = EXP( MIN(1000., rhl_clr) / 100. ) * pdf_e3
         
         IF ( ok_warm_cloud ) THEN
@@ -790,5 +786,5 @@
             dtime, pplay(i), temp(i), qsat(i), qsatl(i), gamma_cond(i), &
             rcont_seri(i), flight_dist(i), cldfra(i), qvc(i), &
-            V_cell, M_cell, pdf_loc, pdf_scale, pdf_alpha, &
+            dz, V_cell, pdf_loc, pdf_scale, pdf_alpha, &
             Tcritcont(i), qcritcont(i), potcontfraP(i), potcontfraNP(i), contfra(i), &
             dcf_avi(i), dqvc_avi(i), dqi_avi(i) &
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90	(revision 5453)
@@ -236,4 +236,7 @@
   REAL, SAVE, PROTECTED :: linear_contrails_lifetime=10800.  ! [s] timescale of the lifetime of linear contrails
   !$OMP THREADPRIVATE(linear_contrails_lifetime)
+
+  REAL, SAVE, PROTECTED :: initial_width_contrails=200.      ! [m] initial width of the linear contrails formed
+  !$OMP THREADPRIVATE(initial_width_contrails)
   !--End of the parameters for aviation
 
@@ -475,4 +478,5 @@
     CALL getin_p('prop_efficiency_aviation',prop_efficiency_aviation)
     CALL getin_p('linear_contrails_lifetime',linear_contrails_lifetime)
+    CALL getin_p('initial_width_contrails',initial_width_contrails)
 
 
@@ -566,4 +570,5 @@
     WRITE(lunout,*) 'lscp_ini, prop_efficiency_aviation:', prop_efficiency_aviation
     WRITE(lunout,*) 'lscp_ini, linear_contrails_lifetime:', linear_contrails_lifetime
+    WRITE(lunout,*) 'lscp_ini, initial_width_contrails:', initial_width_contrails
 
 
Index: LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5453)
@@ -665,4 +665,6 @@
 
 !-- LSCP - aviation and contrails variables
+      REAL, SAVE, ALLOCATABLE :: d_q_avi(:,:)
+      !$OMP THREADPRIVATE(d_q_avi)
       REAL, SAVE, ALLOCATABLE :: rcont_seri(:,:), d_rcont_dyn(:,:)
       !$OMP THREADPRIVATE(rcont_seri, d_rcont_dyn)
@@ -1215,4 +1217,5 @@
 
 !-- LSCP - aviation and contrails variables
+      ALLOCATE(d_q_avi(klon,klev))
       ALLOCATE(rcont_seri(klon,klev), d_rcont_dyn(klon,klev))
       ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
@@ -1621,5 +1624,5 @@
 
 !-- LSCP - aviation and contrails variables
-      DEALLOCATE(rcont_seri, d_rcont_dyn, flight_dist, flight_h2o)
+      DEALLOCATE(d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o)
       DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
       DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
Index: LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5453)
@@ -2168,4 +2168,10 @@
 
 !-- LSCP - aviation variables
+  TYPE(ctrl_out), SAVE :: o_dqavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'dqavi', 'Water vapor emissions from aviation tendency', 'kg/kg/s', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_rcontseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'rcontseri', 'Contrails fraction to total cloud fraction ratio', '-', (/ ('',i=1,10) /))
+  TYPE(ctrl_out), SAVE :: o_drcontdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
+    'drcontdyn', 'Dynamics contrails fraction to total cloud fraction ratio tendency', 's-1', (/ ('',i=1,10) /))
   TYPE(ctrl_out), SAVE :: o_Tcritcont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'Tcritcont', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
Index: LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90	(revision 5453)
@@ -226,4 +226,5 @@
          o_dqvcadj, o_dqvcsub, o_dqvccon, o_dqvcmix, o_qsatl, o_qsati, &
 !-- LSCP - aviation variables
+         o_rcontseri, o_drcontdyn, o_dqavi, &
          o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
          o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
@@ -349,4 +350,5 @@
          dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
          qsatliq, qsatice, &
+         rcont_seri, d_rcont_dyn, d_q_avi, &
          Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
          dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
@@ -2123,4 +2125,6 @@
 !-- LSCP - aviation variables
        IF (ok_plane_contrail) THEN
+         CALL histwrite_phy(o_rcontseri, rcont_seri)
+         CALL histwrite_phy(o_drcontdyn, d_rcont_dyn)
          CALL histwrite_phy(o_flight_dist, flight_dist)
          CALL histwrite_phy(o_Tcritcont, Tcritcont)
@@ -2134,4 +2138,5 @@
        IF (ok_plane_h2o) THEN
          CALL histwrite_phy(o_flight_h2o, flight_h2o)
+         CALL histwrite_phy(o_dqavi, d_q_avi)
        ENDIF
        
Index: LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5452)
+++ LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5453)
@@ -75,4 +75,5 @@
     USE write_field_phy
     use wxios_mod, ONLY: g_ctx, wxios_set_context
+    USE lmdz_aviation, ONLY : read_aviation_emissions, aviation_water_emissions
     USE lmdz_lscp, ONLY : lscp
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
@@ -327,5 +328,5 @@
        dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
        !-- LSCP - aviation and contrails variables
-       rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
+       d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
        Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
        dcf_avi, dqi_avi, dqvc_avi, &
@@ -3863,4 +3864,19 @@
          ratqs,ratqsc,ratqs_inter_,sigma_qtherm)
 
+    !--Read the aviation emissions
+    IF ( ok_plane_h2o .OR. ok_plane_contrail ) THEN
+      CALL read_aviation_emissions(klon, klev, latitude_deg, longitude_deg, pplay, &
+                                   flight_dist, flight_h2o)
+    ENDIF
+
+    !--Add the water emissions from aviation
+    IF ( ok_plane_h2o ) THEN
+       CALL aviation_water_emissions(klon, klev, phys_tstep, paprs, pplay, &
+            t_seri, q_seri, cell_area, flight_h2o, d_q_avi)
+       CALL add_phys_tend(du0, dv0, dt0, d_q_avi, dql0, dqi0, dqbs0, paprs, &
+            'avi', abortphy, flag_inhib_tend, itap, 0)
+       d_q_avi = d_q_avi / phys_tstep
+    ENDIF
+
     !
     ! Appeler le processus de condensation a grande echelle
@@ -3875,5 +3891,4 @@
 
     IF (ok_new_lscp) THEN
-
  
     DO k = 1, klev
@@ -3883,10 +3898,4 @@
       ENDDO
     ENDDO
-
-
-    !--mise à jour de flight_m et flight_h2o dans leur module
-    !IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
-    !  CALL airplane(debut,pphis,pplay,paprs,t_seri)
-    !ENDIF
 
     CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
