Index: LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml
===================================================================
--- LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml	(revision 5597)
+++ LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml	(revision 5598)
@@ -926,5 +926,6 @@
         <field id="potcontfraP"  long_name="Fraction with pontential persistent contrail"    unit="-" />
         <field id="potcontfraNP" long_name="Fraction with potential non-persistent contrail"    unit="-" />
-        <field id="contfra"    long_name="Linear contrail fraction"    unit="kg/kg" />
+        <field id="contfra"    long_name="Linear contrail fraction"    unit="-" />
+        <field id="qice_cont"  long_name="Contrails ice specific humidity seen by radiation"    unit="kg/kg" />
         <field id="dcontfracir"  long_name="Linear contrail fraction to cirrus cloud fraction tendency"    unit="-" />
         <field id="dcfavi"     long_name="Aviation cloud fraction tendency"    unit="s-1" />
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5597)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5598)
@@ -12,6 +12,6 @@
     icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
     !--AB contrails
-    rcontrail, pclc_nocont, pcltau_nocont, pclemi_nocont, pch_nocont, &
-    pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
+    contfra, radocond_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
+    pch_nocont, pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
 
   ! Interface between the LMDZ physics monitor and the cloud properties calculation routines
@@ -94,5 +94,6 @@
 
   !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
-  REAL, INTENT(IN) :: rcontrail(klon, klev)      ! ratio of contrails to total cloud fraction [-]
+  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
+  REAL, INTENT(IN)  :: radocond_cont(klon, klev) ! contrails condensed water seen by radiation [kg/kg]
   REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
   REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
@@ -133,6 +134,6 @@
     icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
     !--AB for contrails
-    rcontrail, pclc_nocont, pcltau_nocont, pclemi_nocont, pch_nocont, &
-    pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
+    contfra, radocond_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
+    pch_nocont, pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
   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 5597)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90	(revision 5598)
@@ -11,6 +11,6 @@
     icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, &
     !--AB contrails
-    rcontrail, pclc_nocont, pcltau_nocont, pclemi_nocont, pch_nocont, &
-    pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
+    contfra, radocond_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
+    pch_nocont, pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
 
   USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc
@@ -118,5 +118,6 @@
 
   !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
-  REAL, INTENT(IN) :: rcontrail(klon, klev)      ! ratio of contrails to total cloud fraction [-]
+  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
+  REAL, INTENT(IN)  :: radocond_cont(klon, klev) ! contrails condensed water seen by radiation [kg/kg]
   REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
   REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
@@ -247,6 +248,6 @@
       DO k = 1, klev
         DO i = 1, klon
-          pclc_nocont(i,k) = pclc(i,k) * ( 1. - rcontrail(i,k) )
-          xfiwc_nocont(i, k) = icefrac_optics(i, k)*radocond(i, k)*( 1. - rcontrail(i,k) )
+          pclc_nocont(i,k) = pclc(i, k) - contfra(i, k)
+          xfiwc_nocont(i, k) = icefrac_optics(i, k) * radocond(i, k) - radocond_cont(i, k)
         ENDDO
       ENDDO
@@ -393,6 +394,6 @@
             !--vs contrail cirrus in the gridbox
             !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
-            rei = rei * ( 1. - rcontrail(i,k) ) &
-                + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
+            rei = ( rei * pclc_nocont(i,k) &
+                + re_ice_crystals_contrails * 1.E6 * contfra(i,k) ) / pclc(i,k)
           ENDIF
           pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
@@ -539,4 +540,5 @@
             reice_nocont(i,k) = 0.
             pclc_nocont(i,k) = 0.
+            pclc(i,k) = contfra(i,k)
             pcltau_nocont(i,k) = 0.
             pclemi_nocont(i, k) = 0.
@@ -552,6 +554,6 @@
           !--vs contrail cirrus in the gridbox
           !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
-          rei = rei * ( 1. - rcontrail(i,k) ) &
-              + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
+          rei = ( rei * pclc_nocont(i,k) &
+              + re_ice_crystals_contrails * 1.E6 * contfra(i,k) ) / pclc(i,k)
         ENDIF
 
@@ -703,11 +705,11 @@
     DO k = klev, 1, -1
       DO i = 1, klon
-        zclear(i) = zclear(i)*(1.-max(pclc(i,k)*rcontrail(i,k),zcloud(i)))/(1.-min(real( &
+        zclear(i) = zclear(i)*(1.-max(contfra(i,k),zcloud(i)))/(1.-min(real( &
           zcloud(i),kind=8),1.-zepsec))
         pct_cont(i) = 1. - zclear(i)
-        zcloud(i) = pclc(i, k)*rcontrail(i,k)
+        zcloud(i) = contfra(i,k)
         IF (paprs(i,k)<prmhc) THEN
-          pch_nocont(i) = pch_nocont(i)*(1.-max(pclc_nocont(i,k),zcloudh(i)))/(1.-min(real(zcloudh &
-            (i),kind=8),1.-zepsec))
+          pch_nocont(i) = pch_nocont(i)*(1.-max(pclc_nocont(i,k),zcloudh(i)))/(1.-min(real( &
+            zcloudh(i),kind=8),1.-zepsec))
           zcloudh(i) = pclc_nocont(i, k)
         ENDIF
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5597)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5598)
@@ -25,5 +25,5 @@
      dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
      rcont_seri, flight_dist, flight_h2o,               &
-     contfra, Tcritcont, qcritcont,                     &
+     contfra, radocond_cont, Tcritcont, qcritcont,      &
      potcontfraP, potcontfraNP, dcontfra_cir, dcf_avi,  &
      dqi_avi, dqvc_avi, cloudth_sth,cloudth_senv,       &
@@ -239,4 +239,5 @@
 
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: contfra        !--linear contrail fraction [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond_cont  !--condensed water in linear contrails used in the radiation scheme [kg/kg]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
@@ -432,4 +433,5 @@
 dqvc_mix(:,:)   = 0.
 contfra(:,:)    = 0.
+radocond_cont(:,:)= 0.
 Tcritcont(:,:)  = missing_val
 qcritcont(:,:)  = missing_val
@@ -1139,4 +1141,5 @@
           !--This barrier should never be activated
           rcont_seri(i,k) = MIN(MAX(rcont_seri(i,k), 0.), 1.)
+          radocond_cont(i,k) = radocond(i,k) * rcont_seri(i,k)
         ENDIF
 
Index: LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5597)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5598)
@@ -679,6 +679,6 @@
       REAL, SAVE, ALLOCATABLE :: potcontfraP(:,:), potcontfraNP(:,:)
       !$OMP THREADPRIVATE(potcontfraP, potcontfraNP)
-      REAL, SAVE, ALLOCATABLE :: contfra(:,:), dcontfra_cir(:,:)
-      !$OMP THREADPRIVATE(contfra, dcontfra_cir)
+      REAL, SAVE, ALLOCATABLE :: contfra(:,:), radocond_cont(:,:), dcontfra_cir(:,:)
+      !$OMP THREADPRIVATE(contfra, radocond_cont, dcontfra_cir)
       REAL, SAVE, ALLOCATABLE :: dcf_avi(:,:), dqi_avi(:,:), dqvc_avi(:,:)
       !$OMP THREADPRIVATE(dcf_avi, dqi_avi, dqvc_avi)
@@ -1246,5 +1246,5 @@
       ALLOCATE(Tcritcont(klon,klev), qcritcont(klon,klev))
       ALLOCATE(potcontfraP(klon,klev), potcontfraNP(klon,klev))
-      ALLOCATE(contfra(klon,klev), dcontfra_cir(klon,klev))
+      ALLOCATE(contfra(klon,klev), radocond_cont(klon,klev), dcontfra_cir(klon,klev))
       ALLOCATE(dcf_avi(klon,klev), dqi_avi(klon,klev), dqvc_avi(klon,klev))
       ALLOCATE(cldfra_nocont(klon,klev), cldtau_nocont(klon,klev), cldemi_nocont(klon,klev))
@@ -1659,5 +1659,5 @@
       DEALLOCATE(d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o)
       DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
-      DEALLOCATE(contfra, dcontfra_cir)
+      DEALLOCATE(contfra, radocond_cont, dcontfra_cir)
       DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
       DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
Index: LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5597)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5598)
@@ -2202,4 +2202,6 @@
   TYPE(ctrl_out), SAVE :: o_contfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'contfra', 'Linear contrail fraction', '-', (/ ('', i=1,10)/))
+  TYPE(ctrl_out), SAVE :: o_qice_cont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qice_cont', 'Linear contrails ice specific humidity seen by radiation', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_dcontfracir = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     'dcontfracir', 'Linear contrail fraction to cirrus cloud fraction tendency', '-', (/ ('', 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 5597)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90	(revision 5598)
@@ -229,5 +229,5 @@
          o_issrfra250to300, o_issrfra300to400, o_issrfra400to500, &
 !-- LSCP - aviation variables
-         o_rcontseri, o_drcontdyn, o_dqavi, o_contfra, &
+         o_rcontseri, o_drcontdyn, o_dqavi, o_contfra, o_qice_cont, &
          o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
          o_dcontfracir, o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
@@ -358,5 +358,5 @@
          issrfra100to150, issrfra150to200, issrfra200to250, &
          issrfra250to300, issrfra300to400, issrfra400to500, &
-         rcont_seri, d_rcont_dyn, d_q_avi, contfra, &
+         rcont_seri, d_rcont_dyn, d_q_avi, contfra, radocond_cont, &
          Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
          dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
@@ -2156,4 +2156,5 @@
          CALL histwrite_phy(o_potcontfraNP, potcontfraNP)
          CALL histwrite_phy(o_contfra, contfra)
+         CALL histwrite_phy(o_qice_cont, radocond_cont)
          CALL histwrite_phy(o_dcontfracir, dcontfra_cir)
          CALL histwrite_phy(o_dcfavi, dcf_avi)
Index: LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5597)
+++ LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5598)
@@ -329,5 +329,5 @@
        !-- LSCP - aviation and contrails variables
        d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
-       contfra, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
+       contfra, radocond_cont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
        dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, &
        cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
@@ -3930,5 +3930,5 @@
          dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
          rcont_seri, flight_dist, flight_h2o, &
-         contfra, Tcritcont, qcritcont, potcontfraP, &
+         contfra, radocond_cont, Tcritcont, qcritcont, potcontfraP, &
          potcontfraNP, dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, &
          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
@@ -4483,6 +4483,6 @@
                zfice, dNovrN, ptconv, rnebcon, clwcon, &
                !--AB contrails
-               rcont_seri, cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
-               contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)
+               contfra, radocond_cont, cldfra_nocont, cldtau_nocont, cldemi_nocont, &
+               cldh_nocont, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)
 
        !