Index: /LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90
===================================================================
--- /LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5487)
+++ /LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5488)
@@ -678,158 +678,4 @@
 !  INCLUDE 'netcdf.inc'
 !
-!  !--------------------------------------------------------
-!  !--input variables 
-!  !--------------------------------------------------------
-!  LOGICAL, INTENT(IN) :: debut
-!  REAL, INTENT(IN)    :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev)
-!
-!  !--------------------------------------------------------
-!  !	... Local variables
-!  !--------------------------------------------------------
-!
-!  CHARACTER (LEN=20) :: modname='airplane_mod'
-!  INTEGER :: i, k, kori, iret, varid, error, ncida, klona
-!  INTEGER,SAVE :: nleva, ntimea
-!!$OMP THREADPRIVATE(nleva,ntimea)
-!  REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:)    !--km/s
-!  REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:)   !--molec H2O/cm3/s
-!  REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:)
-!  REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:)
-!  REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:)
-!!$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta)
-!  REAL :: zalt(klon,klev+1)
-!  REAL :: zrho, zdz(klon,klev), zfrac
-!
-!  !
-!  IF (debut) THEN
-!  !--------------------------------------------------------------------------------
-!  !       ... Open the file and read airplane emissions
-!  !--------------------------------------------------------------------------------
-!  !
-!  IF (is_mpi_root .AND. is_omp_root) THEN
-!      !
-!      iret = nf_open('aircraft_phy.nc', 0, ncida)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1)
-!      ! ... Get lengths
-!      iret = nf_inq_dimid(ncida, 'time', varid)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1)
-!      iret = nf_inq_dimlen(ncida, varid, ntimea)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1)
-!      iret = nf_inq_dimid(ncida, 'vector', varid)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1)
-!      iret = nf_inq_dimlen(ncida, varid, klona)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1)
-!      iret = nf_inq_dimid(ncida, 'lev', varid)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
-!      iret = nf_inq_dimlen(ncida, varid, nleva)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1)
-!      !
-!      IF ( klona /= klon_glo ) THEN
-!        WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo
-!        CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1)
-!      ENDIF
-!      !
-!      IF ( ntimea /= 12 ) THEN
-!        WRITE(lunout,*) 'ntimea=', ntimea
-!        CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1)
-!      ENDIF
-!      !
-!      ALLOCATE(zmida(nleva), STAT=error)
-!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1)
-!      ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error)
-!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1)
-!      ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error)
-!      IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1)
-!      !
-!      iret = nf_inq_varid(ncida, 'lev', varid)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1)
-!      iret = nf_get_var_double(ncida, varid, zmida)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1)
-!      !
-!      iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid)  !--CO2 as a proxy for m flown - 
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1)
-!      iret = nf_get_var_double(ncida, varid, pkm_airpl_glo)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1)
-!      !
-!      iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1)
-!      iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo)
-!      IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1)
-!      ! 
-!     ENDIF    !--is_mpi_root and is_omp_root
-!     !
-!     CALL bcast(nleva)
-!     CALL bcast(ntimea)
-!     ! 
-!     IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error)
-!     IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error)
-!     !
-!     ALLOCATE(pkm_airpl(klon,nleva,ntimea))
-!     ALLOCATE(ph2o_airpl(klon,nleva,ntimea))
-!     !
-!     ALLOCATE(flight_m(klon,klev))
-!     ALLOCATE(flight_h2o(klon,klev))
-!     !
-!     CALL bcast(zmida)
-!     zinta(1)=0.0                                   !--surface
-!     DO k=2, nleva
-!       zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0  !--conversion from km to m
-!     ENDDO
-!     zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface
-!     !print *,'zinta=', zinta
-!     !
-!     CALL scatter(pkm_airpl_glo,pkm_airpl)
-!     CALL scatter(ph2o_airpl_glo,ph2o_airpl)
-!     !
-!!$OMP MASTER
-!     IF (is_mpi_root .AND. is_omp_root) THEN
-!        DEALLOCATE(pkm_airpl_glo)
-!        DEALLOCATE(ph2o_airpl_glo)
-!     ENDIF   !--is_mpi_root
-!!$OMP END MASTER
-!
-!  ENDIF !--debut
-!!
-!!--compute altitude of model level interfaces
-!!
-!  DO i = 1, klon
-!    zalt(i,1)=pphis(i)/RG         !--in m
-!  ENDDO
-!!
-!  DO k=1, klev
-!    DO i = 1, klon
-!      zrho=pplay(i,k)/t_seri(i,k)/RD
-!      zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG
-!      zalt(i,k+1)=zalt(i,k)+zdz(i,k)   !--in m
-!    ENDDO
-!  ENDDO
-!!
-!!--vertical reprojection 
-!!
-!  flight_m(:,:)=0.0 
-!  flight_h2o(:,:)=0.0 
-!!
-!  DO k=1, klev
-!    DO kori=1, nleva 
-!      DO i=1, klon
-!        !--fraction of layer kori included in layer k
-!        zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori))
-!        !--reproject
-!        flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac 
-!        !--reproject
-!        flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac   
-!      ENDDO
-!    ENDDO
-!  ENDDO
-!!
-!  DO k=1, klev
-!    DO i=1, klon
-!      !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell
-!      flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3
-!      flight_m(i,k)=flight_m(i,k)*100.0  !--x100 to augment signal to noise 
-!      !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell
-!      flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6
-!    ENDDO
-!  ENDDO
 !!
 !END SUBROUTINE airplane
Index: /LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90
===================================================================
--- /LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5487)
+++ /LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5488)
@@ -10,5 +10,5 @@
     reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
     reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
-    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon)
+    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail)
 
   ! Interface between the LMDZ physics monitor and the cloud properties calculation routines
@@ -47,4 +47,6 @@
   LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection
   LOGICAL, INTENT(IN) :: ok_newmicro, ok_aie
+
+  REAL, INTENT(IN) :: rcontrail(klon, klev) ! ratio of contrail to total cloud fraction, used only if ok_plane_contrail=y [-]
 
   ! inout:
@@ -116,5 +118,5 @@
     reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
     reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
-    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon)
+    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail)
   ELSE
     CALL nuage (paprs, pplay, &
Index: /LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90
===================================================================
--- /LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90	(revision 5487)
+++ /LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90	(revision 5488)
@@ -9,5 +9,5 @@
     reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
     reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
-    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon)
+    icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, rcontrail)
 
   USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc
@@ -28,4 +28,5 @@
   USE lmdz_cloud_optics_prop_ini , ONLY : ok_icefra_lscp, rei_max, rei_min
   USE lmdz_cloud_optics_prop_ini , ONLY : zepsec, novlp, iflag_ice_thermo, ok_new_lscp
+  USE lmdz_cloud_optics_prop_ini , ONLY : ok_plane_contrail, re_ice_crystals_contrails
   
 
@@ -69,4 +70,6 @@
 
   LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection
+
+  REAL, INTENT(IN) :: rcontrail(klon, klev) ! ratio of contrails to total cloud fraction, used only if ok_plane_contrail=y [-]
 
   ! inout:
@@ -333,4 +336,10 @@
 
           IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
+          IF ( ok_plane_contrail ) THEN
+            !--If contrails are activated, rei is a weighted average between the natural
+            !--rei and the contrails rei, with the weights being the fraction of natural
+            !--vs contrail cirrus in the gridbox
+            rei = rei * ( 1. - rcontrail(i,k) ) + re_ice_crystals_contrails * rcontrail(i,k)
+          ENDIF
           pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
             zfiwp_var*(3.448E-03+2.431/rei)
@@ -442,4 +451,10 @@
         IF (zflwp_var==0.) rel = 1.
         IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
+        IF ( ok_plane_contrail ) THEN
+          !--If contrails are activated, rei is a weighted average between the natural
+          !--rei and the contrails rei, with the weights being the fraction of natural
+          !--vs contrail cirrus in the gridbox
+          rei = rei * ( 1. - rcontrail(i,k) ) + re_ice_crystals_contrails * rcontrail(i,k)
+        ENDIF
         pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
           rei)
Index: /LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop_ini.f90
===================================================================
--- /LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop_ini.f90	(revision 5487)
+++ /LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop_ini.f90	(revision 5488)
@@ -11,4 +11,5 @@
   LOGICAL, PROTECTED :: ok_cdnc
   LOGICAL, PROTECTED :: ok_icefra_lscp, ok_new_lscp
+  LOGICAL, PROTECTED :: ok_plane_contrail
   REAL, PROTECTED :: bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
   REAL, ALLOCATABLE :: latitude_deg(:)
@@ -22,4 +23,5 @@
   REAL, PROTECTED :: rei_max, rei_min
   REAL, PROTECTED :: zepsec
+  REAL, PROTECTED :: re_ice_crystals_contrails=7.5E-6 ! [m] effective radius of ice crystals in contrails
   REAL, PARAMETER :: thres_tau=0.3, thres_neb=0.001
   REAL, PARAMETER :: prmhc=440.*100. ! Pressure between medium and high level cloud in Pa
@@ -39,4 +41,5 @@
 !$OMP THREADPRIVATE(rad_chau1, rad_chau2, rad_froid, rei_max, rei_min)
 !$OMP THREADPRIVATE(zepsec) 
+!$OMP THREADPRIVATE(re_ice_crystals_contrails) 
 
   
@@ -46,5 +49,6 @@
        & ok_cdnc_in, bl95_b0_in, &
        & bl95_b1_in, latitude_deg_in, rpi_in, rg_in, rd_in, zepsec_in, novlp_in, &
-       & iflag_ice_thermo_in, ok_new_lscp_in)
+       & iflag_ice_thermo_in, ok_new_lscp_in, &
+       & ok_plane_contrail_in)
 
     USE ioipsl_getin_p_mod, ONLY : getin_p
@@ -56,4 +60,5 @@
     INTEGER, INTENT(IN) :: novlp_in, iflag_ice_thermo_in
     LOGICAL, INTENT(IN) :: ok_cdnc_in, ok_new_lscp_in
+    LOGICAL, INTENT(IN) :: ok_plane_contrail_in
     REAL, INTENT(IN) :: bl95_b0_in, bl95_b1_in
     REAL, INTENT(IN) :: latitude_deg_in(klon)
@@ -77,4 +82,5 @@
     iflag_ice_thermo = iflag_ice_thermo_in
     ok_new_lscp = ok_new_lscp_in
+    ok_plane_contrail = ok_plane_contrail_in
     
     call getin_p('cdnc_min',cdnc_min)
@@ -98,4 +104,5 @@
     rei_max = 61.29
     CALL getin_p('rei_max',rei_max)
+    CALL getin_p('re_ice_crystals_contrails', re_ice_crystals_contrails)
 
    
Index: /LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5487)
+++ /LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5488)
@@ -1877,5 +1877,6 @@
                                   & ok_cdnc, bl95_b0, &
                                   & bl95_b1, latitude_deg, rpi, rg, rd, &
-                                  & zepsec, novlp, iflag_ice_thermo, ok_new_lscp)
+                                  & zepsec, novlp, iflag_ice_thermo, ok_new_lscp, &
+                                  & ok_plane_contrail)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -4465,5 +4466,5 @@
                ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
                reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
-               zfice, dNovrN, ptconv, rnebcon, clwcon)
+               zfice, dNovrN, ptconv, rnebcon, clwcon, rcont_seri)
 
        !
