Index: LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml
===================================================================
--- LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml	(revision 5585)
+++ LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml	(revision 5589)
@@ -912,8 +912,8 @@
         <field id="qsati"      long_name="Saturation with respect to ice"    unit="kg/kg" />
 
-	<field id="rcontseri"  long_name="Linear contrail fraction to cloud fraction ratio" unit="-" />
-	<field id="drcontdyn"  long_name="Dynamics linear contrail fraction to cloud fraction ratio tendency" unit="s-1" />
-	<field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
-	<field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
+	    <field id="rcontseri"  long_name="Linear contrail fraction to cloud fraction ratio" unit="-" />
+	    <field id="drcontdyn"  long_name="Dynamics linear contrail fraction to cloud fraction ratio tendency" unit="s-1" />
+	    <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
+	    <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
         <field id="potcontfraP"  long_name="Fraction with pontential persistent contrail"    unit="-" />
         <field id="potcontfraNP" long_name="Fraction with potential non-persistent contrail"    unit="-" />
@@ -925,5 +925,18 @@
         <field id="flightdist" long_name="Aviation flown distance concentration"    unit="m/s/m3" />
         <field id="flighth2o"  long_name="Aviation emitted H2O concentration"    unit="kg H2O/s/m3" />
-	<field id="dqavi"      long_name="Water vapor emissions from aviation tendency"    unit="kg/kg/s" />
+	    <field id="dqavi"      long_name="Water vapor emissions from aviation tendency"    unit="kg/kg/s" />
+	    <field id="cldfra_nocont" long_name="Cloud fraction w/o contrails"    unit="-" />
+	    <field id="cldtau_nocont" long_name="Cloud optical thickness w/o contrails"    unit="1" />
+	    <field id="cldemi_nocont" long_name="Cloud optical emissivity w/o contrails"    unit="1" />
+	    <field id="cldh_nocont"   long_name="High-level cloudiness w/o contrails"    unit="-" />
+	    <field id="contcov"       long_name="Total contrails cover"    unit="-" />
+        <field id="iwp_nocont"    long_name="Cloud ice water path w/o contrails"    unit="kg/m2" />
+        <field id="iwc_nocont"    long_name="Cloud ice water content seen by radiation w/o contrails"    unit="kg/m3" />
+	    <field id="ref_ice_nocont" long_name="Effective radius of ice crystals w/o contrails"    unit="microns" />
+        <field id="tops_nocont"   long_name="Solar rad. at TOA w/o contrails"    unit="W/m2" />
+        <field id="topl_nocont"   long_name="IR rad. at TOA w/o contrails"    unit="W/m2" />
+        <field id="sols_nocont"   long_name="Solar rad. at surf. w/o contrails"    unit="W/m2" />
+        <field id="soll_nocont"   long_name="IR rad. at surface w/o contrails"    unit="W/m2" />
+        <field id="nettop_nocont" long_name="Net dn radiatif flux at TOA w/o contrails"    unit="W/m2" />
 
 	<field id="fluxt"     long_name="flux h"     unit="W/m2" />
Index: LMDZ6/branches/contrails/libf/phylmd/clesphys_mod_h.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/clesphys_mod_h.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/clesphys_mod_h.f90	(revision 5589)
@@ -45,5 +45,5 @@
           , iflag_rrtm, ok_strato, ok_hines, ok_qch4                    &
           , iflag_ice_thermo, ok_ice_supersat, ok_no_issr_strato       &
-          , ok_plane_h2o, ok_plane_contrail                            &
+          , ok_plane_h2o, ok_plane_contrail, ok_rad_contrail           &
           , ok_gwd_rando, NSW, iflag_albedo                            &
           , ok_chlorophyll, ok_conserv_q, adjust_tropopause             &
@@ -142,5 +142,5 @@
   INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
   LOGICAL :: ok_ice_supersat, ok_no_issr_strato
-  LOGICAL :: ok_plane_h2o, ok_plane_contrail
+  LOGICAL :: ok_plane_h2o, ok_plane_contrail, ok_rad_contrail
   LOGICAL :: ok_chlorophyll
   LOGICAL :: ok_strato
@@ -203,5 +203,5 @@
   !$OMP      , iflag_rrtm, ok_strato, ok_hines, ok_qch4                    &
   !$OMP      , iflag_ice_thermo, ok_ice_supersat, ok_no_issr_strato       &
-  !$OMP      , ok_plane_h2o, ok_plane_contrail                            &
+  !$OMP      , ok_plane_h2o, ok_plane_contrail, ok_rad_contrail           &
   !$OMP      , ok_gwd_rando, NSW, iflag_albedo                            &
   !$OMP      , ok_chlorophyll, ok_conserv_q, adjust_tropopause             &
Index: LMDZ6/branches/contrails/libf/phylmd/conf_phys_m.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/conf_phys_m.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/conf_phys_m.f90	(revision 5589)
@@ -173,5 +173,5 @@
     INTEGER,SAVE :: iflag_ice_thermo_omp
     LOGICAL,SAVE :: ok_ice_supersat_omp, ok_no_issr_strato_omp
-    LOGICAL,SAVE :: ok_plane_h2o_omp, ok_plane_contrail_omp
+    LOGICAL,SAVE :: ok_plane_h2o_omp, ok_plane_contrail_omp, ok_rad_contrail_omp
     REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
     INTEGER,SAVE :: iflag_sic_omp, iflag_inertie_omp
@@ -1378,4 +1378,12 @@
     ok_plane_contrail_omp = .FALSE.
     CALL getin('ok_plane_contrail',ok_plane_contrail_omp)
+
+    !Config Key  = ok_rad_contrail
+    !Config Desc = double call to radiative routine w/ and w/o contrails
+    !Config Def  = .FALSE.
+    !Config Help =
+    !
+    ok_rad_contrail_omp = .FALSE.
+    CALL getin('ok_rad_contrail',ok_rad_contrail_omp)
 
     !
@@ -2357,4 +2365,5 @@
     ok_plane_h2o = ok_plane_h2o_omp
     ok_plane_contrail = ok_plane_contrail_omp
+    ok_rad_contrail = ok_rad_contrail_omp
     top_height = top_height_omp
     overlap = overlap_omp
@@ -2783,4 +2792,5 @@
     WRITE(lunout,*) ' ok_plane_h2o = ',ok_plane_h2o
     WRITE(lunout,*) ' ok_plane_contrail = ',ok_plane_contrail
+    WRITE(lunout,*) ' ok_rad_contrail = ',ok_rad_contrail
     WRITE(lunout,*) ' overlap = ',overlap 
     WRITE(lunout,*) ' cdmmax = ',cdmmax 
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90	(revision 5589)
@@ -635,6 +635,8 @@
     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 iophy, ONLY : io_lon, io_lat
     USE lmdz_xios
     USE print_control_mod, ONLY: lunout
+    USE readaerosol_mod, ONLY: check_err
     USE lmdz_lscp_ini, ONLY: EI_H2O_aviation
     IMPLICIT NONE
@@ -648,4 +650,21 @@
     INTEGER :: ierr
 
+    ! FOR REGULAR LON LAT
+    CHARACTER(len=30)     :: fname
+    INTEGER               :: nid, dimid, varid
+    INTEGER               :: imth, i, j, k
+    REAL                  :: npole, spole
+    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D
+    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_h2o_glo2D
+    REAL, ALLOCATABLE, DIMENSION(:)       :: vartmp_lev
+    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: vartmp
+    REAL, ALLOCATABLE, DIMENSION(:)       :: lon_src              ! longitudes in file
+    REAL, ALLOCATABLE, DIMENSION(:)       :: lat_src, lat_src_inv ! latitudes in file
+    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
+    INTEGER                               :: nbp_lon, nbp_lat
+
+
+IF (grid_type==unstructured) THEN
+
     ! Get number of vertical levels and level values
     IF (is_omp_master) CALL xios_get_axis_attr( "aviation_lev", n_glo=nleva )
@@ -654,9 +673,9 @@
     ! Allocation of arrays
     ALLOCATE(flight_dist_read(klon, nleva,1), STAT=ierr)
-    IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist',1)
+    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1)
     ALLOCATE(flight_h2o_read(klon, nleva,1), STAT=ierr)
-    IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_h2o',1)
+    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_h2o',1)
     ALLOCATE(aviation_lev(nleva), STAT=ierr)
-    IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate aviation_lev',1)
+    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1)
 
     ! Read the data from the file
@@ -675,8 +694,240 @@
 
     ! 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, flight_dist_read)
     CALL scatter_omp(flight_h2o_mpi, flight_h2o_read)
     CALL bcast_omp(aviation_lev)
+
+ELSE
+
+    nbp_lon=nbp_lon_
+    nbp_lat=nbp_lat_
+    
+    IF (is_mpi_root) THEN
+      ALLOCATE(lon_src(nbp_lon))
+      ALLOCATE(lat_src(nbp_lat))
+      ALLOCATE(lat_src_inv(nbp_lat))
+    ELSE
+      ALLOCATE(flight_dist_glo2D(0,0,0,0))
+      ALLOCATE(flight_h2o_glo2D(0,0,0,0))
+    ENDIF
+
+    IF (is_mpi_root .AND. is_omp_root) THEN
+
+! 1) Open file 
+!****************************************************************************************
+      fname = 'aviation.nc'
+      CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, nid), "pb open "//TRIM(fname) )
+
+
+! Test for equal longitudes and latitudes in file and model
+!****************************************************************************************
+      ! Read and test longitudes
+      CALL check_err( nf90_inq_varid(nid, 'longitude', varid), "pb inq lon" )
+      CALL check_err( nf90_get_var(nid, varid, lon_src(:)), "pb get lon" )
+      
+      IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
+         WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
+         WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
+         WRITE(lunout,*) 'longitudes in model :', io_lon
+       
+         CALL abort_physic('lmdz_aviation', 'longitudes are not the same in file and model',1)
+      END IF
+
+      ! Read and test latitudes
+      CALL check_err( nf90_inq_varid(nid, 'latitude', varid),"pb inq lat" )
+      CALL check_err( nf90_get_var(nid, varid, lat_src(:)),"pb get lat" )
+
+      ! Invert source latitudes
+      DO j = 1, nbp_lat
+         lat_src_inv(j) = lat_src(nbp_lat+1-j)
+      END DO
+
+      IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
+         ! Latitudes are the same
+         invert_lat=.FALSE.
+      ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
+         ! Inverted source latitudes correspond to model latitudes
+         WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
+         invert_lat=.TRUE.
+      ELSE
+         WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
+         WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src      
+         WRITE(lunout,*) 'latitudes in model :', io_lat
+         CALL abort_physic('lmdz_aviation', 'latitudes do not correspond between file and model',1)
+      END IF
+
+       
+! 2) Find vertical dimension nleva
+!****************************************************************************************
+       ierr = nf90_inq_dimid(nid, 'pressure_Pa', dimid) 
+       CALL check_err( nf90_inquire_dimension(nid, dimid, len = nleva), "pb inq dim for PRESNIVS or lev" )
+       
+     ! Allocate variables depending on the number of vertical levels
+       ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_h2o_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1)
+
+       ALLOCATE(aviation_lev(nleva), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 2',1)
+
+! 3) Read all variables from file
+!****************************************************************************************
+
+       ! Get variable id
+       CALL check_err( nf90_inq_varid(nid, 'seg_length_m', varid),"pb inq var seg_length_m" )
+       ! Get the variable
+       CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m"  )
+
+!       ! Get variable id
+!       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
+!       ! Get the variable
+!       CALL check_err( nf90_get_var(nid, varid, flight_h2o_glo2D),"pb get var fuel_burn"  )
+       flight_h2o_glo2D(:,:,:,:) = 0.
+
+       ! Get variable id
+       CALL check_err( nf90_inq_varid(nid, "pressure_Pa", varid),"pb inq var pressure_Pa" )
+       ! Get the variable
+       CALL check_err( nf90_get_var(nid, varid, aviation_lev),"pb get var pressure_Pa" )
+          
+
+! 4) Close file  
+!****************************************************************************************
+       CALL check_err( nf90_close(nid), "pb in close" )
+     
+
+! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
+!****************************************************************************************
+! Test if vertical levels have to be inversed
+
+       IF (aviation_lev(1) < aviation_lev(nleva)) THEN
+          
+          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva), vartmp_lev(nleva))
+
+          ! Inverse vertical levels for flight_dist_glo2D
+          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1)
+          DO k=1, nleva
+             DO j=1, nbp_lat
+                DO i=1, nbp_lon
+                   flight_dist_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
+                END DO
+             END DO
+          END DO
+
+          ! Inverse vertical levels for flight_h2o_glo2D
+          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1)
+          DO k=1, nleva
+             DO j=1, nbp_lat
+                DO i=1, nbp_lon
+                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
+                END DO
+             END DO
+          END DO
+           
+          ALLOCATE(vartmp_lev(nleva))
+          ! Inverte vertical axes for aviation_lev
+          vartmp_lev(:) = aviation_lev(:)
+          DO k=1, nleva
+             aviation_lev(k) = vartmp_lev(nleva+1-k)
+          END DO
+
+          DEALLOCATE(vartmp, vartmp_lev)
+
+       ELSE 
+          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
+       END IF
+
+
+!     - Invert latitudes if necessary
+       IF (invert_lat) THEN
+
+          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva))
+
+          ! Invert latitudes for the variable
+          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) ! use varmth temporarly
+          DO k=1,nleva
+             DO j=1,nbp_lat
+                DO i=1,nbp_lon
+                   flight_dist_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
+                END DO
+             END DO
+          END DO
+
+          ! Invert latitudes for the variable
+          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1) ! use varmth temporarly
+          DO k=1,nleva
+             DO j=1,nbp_lat
+                DO i=1,nbp_lon
+                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
+                END DO
+             END DO
+          END DO
+
+          DEALLOCATE(vartmp)
+       END IF ! invert_lat
+        
+       ! Do zonal mead at poles and distribut at whole first and last latitude
+       DO k=1, nleva
+          npole=0.  ! North pole, j=1
+          spole=0.  ! South pole, j=nbp_lat        
+          DO i=1,nbp_lon
+             npole = npole + flight_dist_glo2D(i,1,k,1)
+             spole = spole + flight_dist_glo2D(i,nbp_lat,k,1)
+          END DO
+          npole = npole/REAL(nbp_lon)
+          spole = spole/REAL(nbp_lon)
+          flight_dist_glo2D(:,1,    k,1) = npole
+          flight_dist_glo2D(:,nbp_lat,k,1) = spole
+       END DO
+
+       ! Do zonal mead at poles and distribut at whole first and last latitude
+       DO k=1, nleva
+          npole=0.  ! North pole, j=1
+          spole=0.  ! South pole, j=nbp_lat        
+          DO i=1,nbp_lon
+             npole = npole + flight_h2o_glo2D(i,1,k,1)
+             spole = spole + flight_h2o_glo2D(i,nbp_lat,k,1)
+          END DO
+          npole = npole/REAL(nbp_lon)
+          spole = spole/REAL(nbp_lon)
+          flight_h2o_glo2D(:,1,    k,1) = npole
+          flight_h2o_glo2D(:,nbp_lat,k,1) = spole
+       END DO
+     
+       ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_h2o_mpi(klon_glo, nleva, 1), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1)
+     
+       ! Transform from 2D to 1D field
+       CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi)
+       CALL grid2Dto1D_glo(flight_h2o_glo2D, flight_h2o_mpi)
+    
+    ELSE
+        ALLOCATE(flight_dist_mpi(0,0,0), flight_h2o_mpi(0,0,0))
+    END IF ! is_mpi_root .AND. is_omp_root
+
+!$OMP BARRIER
+  
+! 6) Distribute to all processes
+!    Scatter global field(klon_glo) to local process domain(klon)
+!    and distribute nleva to all processes
+!****************************************************************************************
+
+    ! Distribute nleva
+    CALL bcast(nleva)
+
+    ! Allocate and distribute pt_ap and pt_b
+    IF (.NOT. ALLOCATED(aviation_lev)) THEN  ! if pt_ap is allocated also pt_b is allocated
+       ALLOCATE(aviation_lev(nleva), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 4',1)
+    END IF
+    CALL bcast(aviation_lev)
+
+    ! Allocate space for output pointer variable at local process
+    ALLOCATE(flight_dist_read(klon, nleva, 1), flight_h2o_read(klon, nleva, 1), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1)
+
+    ! Scatter global field to local domain at local process
+    CALL scatter(flight_dist_mpi, flight_dist_read)
+    CALL scatter(flight_h2o_mpi, flight_h2o_read)
+
+ENDIF
 
 END SUBROUTINE read_aviation_emissions
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90	(revision 5589)
@@ -10,5 +10,8 @@
     reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
     reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
-    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail)
+    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)
 
   ! Interface between the LMDZ physics monitor and the cloud properties calculation routines
@@ -47,6 +50,4 @@
   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:
@@ -92,4 +93,16 @@
   REAL, INTENT(INOUT) :: icefrac_optics(klon, klev)! ice fraction in clouds seen by radiation [-]
 
+  !--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(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
+  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
+  REAL, INTENT(OUT) :: xfiwp_nocont(klon)        ! ice water path (seen by radiation) without contrails [kg/m2]
+  REAL, INTENT(OUT) :: xfiwc_nocont(klon, klev)  ! ice water content seen by radiation without contrails [kg/kg]
+  REAL, INTENT(OUT) :: pclc_nocont(klon, klev)   ! cloud fraction for radiation without contrails [-]
+  REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [m]
+  REAL, INTENT(OUT) :: pclemi_nocont(klon, klev) ! cloud emissivity without contrails [-]
+  REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [m]
+  !--AB
+
   ! Local variables
   !----------------
@@ -118,5 +131,8 @@
     reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
     reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
-    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail)
+    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)
   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 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90	(revision 5589)
@@ -9,5 +9,8 @@
     reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
     reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
-    icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, rcontrail)
+    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)
 
   USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc
@@ -72,6 +75,4 @@
   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:
   REAL, INTENT(INOUT) :: pclc(klon, klev) ! cloud fraction for radiation [-]
@@ -116,4 +117,16 @@
   REAL, INTENT(INOUT) :: icefrac_optics(klon, klev)! ice fraction in clouds seen by radiation [-]
 
+  !--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(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
+  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
+  REAL, INTENT(OUT) :: xfiwp_nocont(klon)        ! ice water path (seen by radiation) without contrails [kg/m2]
+  REAL, INTENT(OUT) :: xfiwc_nocont(klon, klev)  ! ice water content seen by radiation without contrails [kg/kg]
+  REAL, INTENT(OUT) :: pclc_nocont(klon, klev)   ! cloud fraction for radiation without contrails [-]
+  REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [m]
+  REAL, INTENT(OUT) :: pclemi_nocont(klon, klev) ! cloud emissivity without contrails [-]
+  REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [m]
+  !--AB
+
   ! Local variables
   !----------------
@@ -173,4 +186,10 @@
   xflwc = 0.D0
   xfiwc = 0.D0
+  !--AB
+  IF (ok_plane_contrail) THEN
+    xfiwp_nocont = 0.D0
+    xfiwc_nocont = 0.D0
+    reice_nocont = 0.
+  ENDIF
 
   reliq = 0.
@@ -222,4 +241,15 @@
       ENDDO
     ENDDO
+
+    !--AB
+    IF (ok_plane_contrail) THEN
+      !--If contrails are activated, we diagnose iwc without contrails
+      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) )
+        ENDDO
+      ENDDO
+    ENDIF
   ENDIF
 
@@ -363,5 +393,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 * ( 1. - rcontrail(i,k) ) &
+                + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
           ENDIF
           pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
@@ -415,4 +446,12 @@
         pcltau(i, k) = 0.0
         pclemi(i, k) = 0.0
+
+        !--AB
+        IF ( ok_plane_contrail ) THEN
+          reice_nocont(i,k) = 0.
+          pclc_nocont(i,k) = 0.
+          pcltau_nocont(i,k) = 0.
+          pclemi_nocont(i,k) = 0.
+        ENDIF
 
       ELSE
@@ -495,10 +534,26 @@
         IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
         IF ( ok_plane_contrail ) THEN
+          !--Diagnostics of clouds emissivity, optical depth and ice crystal radius
+          !--without contrails
+          IF ( pclc_nocont(i,k) .LE. seuil_neb ) THEN
+            reice_nocont(i,k) = 0.
+            pclc_nocont(i,k) = 0.
+            pcltau_nocont(i,k) = 0.
+            pclemi_nocont(i, k) = 0.
+          ELSE
+            reice_nocont(i,k) = rei
+            pcltau_nocont(i,k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/rei)
+            k_ice = k_ice0 + 1.0/rei
+            pclemi_nocont(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
+          ENDIF
+
           !--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
           !--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)
-        ENDIF
+          rei = rei * ( 1. - rcontrail(i,k) ) &
+              + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
+        ENDIF
+
         pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
           rei)
@@ -519,4 +574,8 @@
       xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
       xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
+      !--AB
+      IF ( ok_plane_contrail ) THEN
+          xfiwp_nocont(i) = xfiwp_nocont(i) + xfiwc_nocont(i, k)*rhodz(i, k)
+      ENDIF
 
     ENDDO
@@ -630,4 +689,34 @@
     pcl(i) = 1. - pcl(i)
   ENDDO
+
+  !--AB we recompute high cloud cover with contrails only
+  !--and without contrails
+  IF ( ok_plane_contrail ) THEN
+    DO i = 1, klon
+      zclear(i) = 1.
+      zcloud(i) = 0.
+      zcloudh(i) = 0.
+      pch_nocont(i) = 1.
+      pct_cont(i) = 1.
+    ENDDO
+
+    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( &
+          zcloud(i),kind=8),1.-zepsec))
+        pct_cont(i) = 1. - zclear(i)
+        zcloud(i) = pclc(i, k)*rcontrail(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))
+          zcloudh(i) = pclc_nocont(i, k)
+        ENDIF
+      ENDDO
+    ENDDO
+
+    DO i = 1, klon
+      pch_nocont(i) = 1. - pch_nocont(i)
+    ENDDO
+  ENDIF ! ok_plane_contrail
 
   ! ========================================================
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90	(revision 5589)
@@ -31,5 +31,5 @@
      qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
      dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
-     dqsmelt, dqsfreez, dcfres, dqsres, dqvcres)
+     dqsmelt, dqsfreez, dqised, dcfsed, dqvcsed)
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -264,7 +264,7 @@
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcfres         !--cloud fraction tendency due to resuspension of ice crystals [kg/kg/s]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsres         !--snow tendency due to resuspension of ice crystals [kg/kg/s]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvcres        !--cloud water vapor tendency due to resuspension of ice crystals [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqised         !--ice water content tendency due to sedmentation of ice crystals [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
 
   ! for thermals
@@ -306,4 +306,5 @@
   REAL, DIMENSION(klon) :: ztupnew
   REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
+  REAL, DIMENSION(klon) :: zqice_sedim ! for sedimentation of ice crystals
   REAL, DIMENSION(klon) :: zrflclr, zrflcld
   REAL, DIMENSION(klon) :: ziflclr, ziflcld
@@ -474,4 +475,5 @@
 zqupnew(:)    = 0.
 zqvapclr(:)   = 0.
+zqice_sedim(:)= 0.
 
 
@@ -529,10 +531,10 @@
       CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
                         zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
-                        zqvapclr, zqupnew, &
+                        zqvapclr, zqupnew, zqice_sedim, &
                         cldfra_in, qvc_in, qliq_in, qice_in, &
                         zrfl, zrflclr, zrflcld, &
                         zifl, ziflclr, ziflcld, &
                         dqreva(:,k), dqssub(:,k), &
-                        dqsres(:,k), dcfres(:,k), dqvcres(:,k) &
+                        dcfsed(:,k), dqvcsed(:,k) &
                         )
 
@@ -974,5 +976,5 @@
                             ctot_vol(:,k), ptconv(:,k), &
                             zt, zq, zoliql, zoliqi, zfice, &
-                            rneb(:,k), znebprecipclr, znebprecipcld, &
+                            rneb(:,k), zqice_sedim, znebprecipclr, znebprecipcld, &
                             zrfl, zrflclr, zrflcld, &
                             zifl, ziflclr, ziflcld, &
@@ -980,5 +982,5 @@
                             dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
                             dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
-                            dqsmelt(:,k), dqsfreez(:,k) &
+                            dqsmelt(:,k), dqsfreez(:,k), dqised(:,k) &
                             )
       DO i = 1, klon
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90	(revision 5589)
@@ -119,9 +119,9 @@
 USE lmdz_lscp_ini,        ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
 USE lmdz_lscp_ini,        ONLY: eps, temp_nowater, ok_weibull_warm_clouds
-USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
+USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds
 USE lmdz_lscp_ini,        ONLY: ok_plane_contrail, ok_no_issr_strato
 
 USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, &
-                                mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
+                                beta_pdf_lscp, temp_thresh_pdf_lscp, &
                                 std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, &
                                 coef_mixing_lscp, coef_shear_lscp, chi_mixing_lscp
@@ -361,4 +361,9 @@
           qvc(i) = 0.
 
+          !--If all the ice has been sublimated, we sublimate
+          !--completely the cloud and do not activate the sublimation
+          !--process
+          qvapincld_new = 0.
+
         !--Else, the cloud is adjusted and sublimated
         ELSE
@@ -406,53 +411,16 @@
           !------------------------------------
 
-          !--If the vapor in cloud is below vapor needed for the cloud to survive
-          IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
-            !--Sublimation of the subsaturated cloud
-            !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
-            !--to use.
-            !--iflag = 1 --> uniform distribution
-            !--iflag = 2 --> exponential distribution
-            !--iflag = 3 --> gamma distribution (Karcher et al 2018)
-
-            IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
-              !--Uniform distribution starting at qvapincld
-              pdf_e1 = 1. / ( 2. * qiceincld )
-
-              dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
-              dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
-                                       * pdf_e1 / 2.
-
-            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
-              !--Exponential distribution starting at qvapincld
-              pdf_alpha = 1. / qiceincld
-              pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
-
-              dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
-              dqt_sub    = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
-                                       + qvapincld - qvapincld_new * pdf_e1 )
-
-            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
-              !--Gamma distribution starting at qvapincld
-              pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
-              pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
-              pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
-              pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
-
-              dcf_sub(i) = - cldfra(i) * pdf_e1
-              dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
-
-            ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN
-              !--Normal distribution around qvapincld + qiceincld with width
-              !--constant in the RHi space
-              pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
-                    / ( std_subl_pdf_lscp / 100. * qsat(i))
-              pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
-              pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
-
-              dcf_sub(i) = - cldfra(i) * pdf_e1
-              dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
-                                         - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
-                                                                                       
-            ENDIF
+          !--If the dissipation process must be activated
+          IF ( qvapincld_new .GT. 0. ) THEN
+            !--Normal distribution around qvapincld + qiceincld with width
+            !--constant in the RHi space
+            pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
+                  / ( std_subl_pdf_lscp / 100. * qsat(i))
+            pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
+            pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
+
+            dcf_sub(i) = - cldfra(i) * pdf_e1
+            dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
+                                       - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
 
             !--Tendencies and diagnostics
@@ -464,5 +432,5 @@
             qvc(i)    = MAX(0., qvc(i) + dqvc_sub(i))
 
-          ENDIF ! qvapincld .LT. qvapincld_new
+          ENDIF ! qvapincld_new .GT. 0.
 
         ENDIF   ! qiceincld .LT. eps
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90	(revision 5589)
@@ -162,7 +162,4 @@
   !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
 
-  INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=4       ! iflag for the distribution of water inside ice clouds
-  !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)
-
   REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
   !$OMP THREADPRIVATE(depo_coef_cirrus)
@@ -173,7 +170,4 @@
   REAL, SAVE, PROTECTED :: std_subl_pdf_lscp=2.              ! [%] standard deviation of the gaussian distribution of water inside ice clouds
   !$OMP THREADPRIVATE(std_subl_pdf_lscp)
-
-  REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] shape factor of the gamma distribution of water inside ice clouds
-  !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
   
   REAL, SAVE, PROTECTED :: beta_pdf_lscp=1.E-3               ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
@@ -260,4 +254,7 @@
   !$OMP THREADPRIVATE(ok_corr_vap_evasub)
 
+  LOGICAL, SAVE, PROTECTED :: ok_growth_precip_deposition=.FALSE.
+  !$OMP THREADPRIVATE(ok_growth_precip_deposition)
+
   REAL, SAVE, PROTECTED :: cld_lc_lsc_snow=2.e-5            ! snow autoconversion coefficient, stratiform. default from  Chaboureau and PInty 2006
   !$OMP THREADPRIVATE(cld_lc_lsc_snow)
@@ -311,4 +308,7 @@
   !$OMP THREADPRIVATE(r_snow)
 
+  REAL, SAVE, PROTECTED :: expo_tau_auto_snow=0.1
+  !$OMP THREADPRIVATE(expo_tau_auto_snow)
+
   REAL, SAVE, PROTECTED :: tau_auto_snow_min=100.           ! Snow autoconversion minimal timescale (when liquid) [s]
   !$OMP THREADPRIVATE(tau_auto_snow_min)
@@ -344,9 +344,9 @@
   !$OMP THREADPRIVATE(snow_fallspeed_cld)
 
-  LOGICAL, SAVE, PROTECTED :: ok_precip_resuspension=.FALSE.! Flag to activate the resuspension of ice crystals
-  !$OMP THREADPRIVATE(ok_precip_resuspension)
-
-  REAL, SAVE, PROTECTED :: snow_thresh_resuspension=1.e-5   ! [kg/m2/s] Threshold flux below which part of the snow flux will be resuspended
-  !$OMP THREADPRIVATE(snow_thresh_resuspension)
+  LOGICAL, SAVE, PROTECTED :: ok_ice_sedim=.FALSE.          ! Flag to activate the sedimentation of ice crystals
+  !$OMP THREADPRIVATE(ok_ice_sedim)
+
+  REAL, SAVE, PROTECTED :: ice_fallspeed_sedim=0.1          ! Ice fallspeed velocity for sedimentation [m/s]
+  !$OMP THREADPRIVATE(ice_fallspeed_sedim)
   !--End of the parameters for poprecip
 
@@ -448,4 +448,5 @@
     CALL getin_p('ok_poprecip',ok_poprecip)
     CALL getin_p('ok_corr_vap_evasub',ok_corr_vap_evasub)
+    CALL getin_p('ok_growth_precip_deposition',ok_growth_precip_deposition)
     CALL getin_p('rain_int_min',rain_int_min) 
     CALL getin_p('gamma_agg',gamma_agg)
@@ -467,14 +468,12 @@
     CALL getin_p('snow_fallspeed_clr',snow_fallspeed_clr)
     CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
-    CALL getin_p('ok_precip_resuspension',ok_precip_resuspension)
-    CALL getin_p('snow_thresh_resuspension',snow_thresh_resuspension)
+    CALL getin_p('ok_ice_sedim',ok_ice_sedim)
+    CALL getin_p('ice_fallspeed_sedim',ice_fallspeed_sedim)
     ! for condensation and ice supersaturation
     CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
     CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds)
-    CALL getin_p('iflag_cloud_sublim_pdf',iflag_cloud_sublim_pdf)
     CALL getin_p('depo_coef_cirrus',depo_coef_cirrus)
     CALL getin_p('capa_cond_cirrus',capa_cond_cirrus)
     CALL getin_p('std_subl_pdf_lscp',std_subl_pdf_lscp)
-    CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
     CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
     CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
@@ -553,4 +552,5 @@
     WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
     WRITE(lunout,*) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
+    WRITE(lunout,*) 'lscp_ini, ok_growth_precip_deposition', ok_growth_precip_deposition
     WRITE(lunout,*) 'lscp_ini, rain_int_min:', rain_int_min
     WRITE(lunout,*) 'lscp_ini, gamma_agg:', gamma_agg
@@ -566,15 +566,13 @@
     WRITE(lunout,*) 'lscp_ini, snow_fallspeed_clr:', snow_fallspeed_clr
     WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
-    WRITE(lunout,*) 'lscp_ini, ok_precip_resuspension:', ok_precip_resuspension
-    WRITE(lunout,*) 'lscp_ini, snow_thresh_resuspension:', snow_thresh_resuspension
+    WRITE(lunout,*) 'lscp_ini, ok_ice_sedim:', ok_ice_sedim
+    WRITE(lunout,*) 'lscp_ini, ice_fallspeed_sedim:', ice_fallspeed_sedim
     ! for condensation and ice supersaturation
     WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
     WRITE(lunout,*) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds
     WRITE(lunout,*) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
-    WRITE(lunout,*) 'lscp_ini, iflag_cloud_sublim_pdf:', iflag_cloud_sublim_pdf
     WRITE(lunout,*) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
     WRITE(lunout,*) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus
     WRITE(lunout,*) 'lscp_ini, std_subl_pdf_lscp:', std_subl_pdf_lscp
-    WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
     WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
     WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp
Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90	(revision 5589)
@@ -540,5 +540,5 @@
           IF (zoliqi(i) .GT. 0.) THEN
             zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & 
-              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
+              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)/zneb(i)*zrho(i)*ffallv)**(-1./dice_velo))
           ELSE
             zfroi=0.
@@ -713,8 +713,8 @@
 SUBROUTINE poprecip_precld( &
            klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
-           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
+           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, qice_sedim, &
            cldfra, qvc, qliq, qice, &
            rain, rainclr, raincld, snow, snowclr, snowcld, &
-           dqreva, dqssub, dqsres, dcfres, dqvcres &
+           dqreva, dqssub, dcfsed, dqvcsed &
            )
 
@@ -722,6 +722,6 @@
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
+USE lmdz_lscp_ini, ONLY : ok_growth_precip_deposition, ok_ice_sedim
 USE lmdz_lscp_ini, ONLY : eps, temp_nowater
-USE lmdz_lscp_ini, ONLY : ok_precip_resuspension, snow_thresh_resuspension
 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
 
@@ -749,4 +749,5 @@
 REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
 REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: qice_sedim     !--ice water content that sedimentated from the layer above [kg/kg]
 
 REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
@@ -765,7 +766,6 @@
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsres         !--snow tendency due to resuspension of ice crystals [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dcfres         !--cloud fraction tendency due to resuspension of ice crystals [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqvcres        !--cloud water vapor tendency due to resuspension of ice crystals [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [s-1]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
 
 
@@ -791,6 +791,4 @@
 !--Specific heat constant
 REAL :: cpair, cpw
-!--For resuspension of ice crystals
-REAL :: dsnowres, qiceinprecip
 
 !--Initialisation
@@ -798,10 +796,7 @@
 dqreva(:) = 0.
 dqssub(:) = 0.
-dqrevap   = 0.
-dqssubl   = 0.
-IF ( ok_precip_resuspension ) THEN
-  dqsres(:) = 0.
-  dcfres(:) = 0.
-  dqvcres(:) = 0.
+IF ( ok_ice_sedim .AND. ok_ice_supersat ) THEN
+  dcfsed(:) = 0.
+  dqvcsed(:) = 0.
 ENDIF
 
@@ -891,4 +886,6 @@
 DO i = 1, klon
 
+  dqrevap = 0.
+  dqssubl = 0.
   !--If there is precipitation from the layer above
   IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
@@ -924,4 +921,11 @@
       IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
         qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
+      ELSE
+        qvapclr = 0.
+      ENDIF
+      IF (  precipfraccld(i)  .GT. eps ) THEN
+        qvapcld = MAX(qtotupnew(i)-qvapclrup(i) , 0.) / qtotupnew(i) * qvap(i) /  precipfraccld(i)
+      ELSE
+        qvapcld = 0.
       ENDIF
     ELSE
@@ -929,5 +933,46 @@
       !--for the evap / subl process
       qvapclr = qvap(i)
-    ENDIF
+      qvapcld = qvap(i)
+    ENDIF
+
+
+    !--------------------
+    !-- ICE SEDIMENTATION
+    !--------------------
+
+    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
+    !--added to the total water content of the gridbox
+    IF ( ok_ice_sedim .AND. ( qice_sedim(i) .GT. 0. ) ) THEN
+      !--Vapor is updated after evaporation/sublimation (it is increased)
+      qvap(i) = qvap(i) + qice_sedim(i)
+      !--Air and precip temperature (i.e., gridbox temperature)
+      !--is updated due to latent heat cooling
+      temp(i) = temp(i) &
+              - qice_sedim(i) * RLSTT / RCPD &
+              / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
+    
+      IF ( ok_ice_supersat ) THEN
+        !--If ice supersaturation is activated, the cloud properties are prognostic.
+        !--The falling ice is then considered a new cloud in the gridbox.
+        !--BEWARE with this parameterisation, we can create a new cloud with a much
+        !--different ice water content and water vapor content than the existing cloud
+        !--(if it exists). This implies than unphysical fluxes of ice and vapor
+        !--occur between the existing cloud and the newly formed cloud (from sedimentation).
+        !--Note also that currently, the precipitation scheme assume a maximum
+        !--random overlap, meaning that the new formed clouds will not be affected
+        !--by vertical wind shear.
+        ! Maybe we should reduce dcfsed?
+        dcfsed(i) = precipfracclr_tmp(i)
+        dqvcsed(i) = qvapclr * dcfsed(i)
+        !--Add tendencies
+        cldfra(i) = cldfra(i) + dcfsed(i)
+        qvc(i) = qvc(i) + dqvcsed(i)
+        qice(i) = qice(i) + qice_sedim(i)
+        !--Normalisation
+        dcfsed(i) = dcfsed(i) / dtime
+        dqvcsed(i) = dqvcsed(i) / dtime
+      ENDIF
+    ENDIF ! ok_ice_sedim
+
 
     !---------------------------
@@ -997,4 +1042,9 @@
       ENDIF
 
+    ELSE
+      !--All the precipitation is sublimated if the fraction is zero
+      drainclreva = - rainclr_tmp(i)
+      dsnowclrsub = - snowclr_tmp(i)
+
     ENDIF ! precipfracclr_tmp .GT. eps
 
@@ -1004,5 +1054,5 @@
     !---------------------------
 
-    IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
+    IF ( ok_growth_precip_deposition .AND. ( temp(i) .LE. RTT ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
       !--Evaporation of liquid precipitation coming from above
       !--in the cloud only
@@ -1011,13 +1061,12 @@
       !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
       !--which does not need a barrier on raincld, because included in the formula
-      !draincldeva = precipfraccld_tmp(i) * MAX(0., &
-      !            - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
-      !            + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
-      !            )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
+      draincldeva = precipfraccld_tmp(i) * MAX(0., &
+                  - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
+                  + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
+                  )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
                
       !--Evaporation is limited by 0
       !--NB. with ok_ice_supersat activated, this barrier should be useless
-      !draincldeva = MIN(0., draincldeva)
-      draincldeva = 0.
+      draincldeva = MIN(0., draincldeva)
 
 
@@ -1081,50 +1130,4 @@
     raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
     snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)
-    
-
-    !--If ok_precip_to_clouds is set to .TRUE., some of the precipitations can go back to clouds
-    !--Only the snow flux can go back to cloud, and only if there is no rain flux
-    !--(mixed-phase clouds)
-    IF ( ok_precip_resuspension ) THEN
-
-      IF ( (snowclr_tmp(i) .GT. 0.) .AND. (precipfracclr_tmp(i) .GT. eps) &
-              .AND. (rainclr_tmp(i) .LT. eps) ) THEN
-        dsnowres = snowclr_tmp(i) * EXP( - snowclr_tmp(i) / precipfracclr_tmp(i) &
-                                           / snow_thresh_resuspension )
-        IF ( dsnowres .GT. 0. ) THEN
-          dqsres(i) = dqsres(i) - dsnowres / dhum_to_dflux(i)
-          !--The following line determines the in-cloud ice water content of the newly formed cloud
-          !--It is a linear combination between the in-precip ice water content (left term) and the
-          !--in-existing-cloud ice water content (right term). This is done so that the existing
-          !--cloud is not too modified, i.e. unphysical ice fluxes going from the existing cloud
-          !--to the new cloud. If the existing cloud is already large, the newly formed cloud
-          !--fraction is reduced to match the in-cloud ice water content of the existing cloud.
-          !--NB. this is done on qi, not on qvc. Not sure what impact that implies
-          qiceinprecip = - dqsres(i) / precipfracclr_tmp(i) * ( 1. - cldfra(i) ) + qice(i)
-          dcfres(i) = - dqsres(i) / qiceinprecip
-          dqvcres(i) = qvapclr * dcfres(i)
-          !--Add tendencies
-          snowcld_tmp(i) = snowcld_tmp(i) + dsnowres
-        ENDIF
-      ENDIF
-
-      IF ( (snowcld_tmp(i) .GT. 0.) .AND. (precipfraccld_tmp(i) .GT. 0.) &
-              .AND. (raincld_tmp(i) .LT. eps) ) THEN
-        dsnowres = snowcld_tmp(i) * EXP( - snowcld_tmp(i) / precipfraccld_tmp(i) &
-                                           / snow_thresh_resuspension )
-        IF ( dsnowres .GT. 0. ) THEN
-          dqsres(i) = dqsres(i) - dsnowres / dhum_to_dflux(i)
-          !--Add tendencies
-          snowcld_tmp(i) = snowcld_tmp(i) + dsnowres
-        ENDIF
-      ENDIF
-
-      !--Add tendencies
-      cldfra(i) = cldfra(i) + dcfres(i)
-      qvc(i) = qvc(i) + dqvcres(i)
-      qice(i) = qice(i) - dqsres(i)
-      !--Recall that qvap is here the total water in the gridbox
-      qvap(i) = qvap(i) - dqsres(i)
-    ENDIF
 
 
@@ -1323,9 +1326,9 @@
 SUBROUTINE poprecip_postcld( &
            klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
-           temp, qvap, qliq, qice, icefrac, cldfra, &
+           temp, qvap, qliq, qice, icefrac, cldfra, qice_sedim, &
            precipfracclr, precipfraccld, &
            rain, rainclr, raincld, snow, snowclr, snowcld, &
            qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
-           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
+           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez, dqised)
 
 USE lmdz_lscp_ini, ONLY : prt_level, lunout
@@ -1338,10 +1341,11 @@
                           rho_rain, r_rain, r_snow, rho_ice,                   &
                           tau_auto_snow_min, tau_auto_snow_max,                &
-                          thresh_precip_frac, eps,                             &
+                          expo_tau_auto_snow, thresh_precip_frac, eps,         &
                           gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
                           iflag_cloudth_vert, iflag_rain_incloud_vol,          &
                           cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
                           rain_fallspeed_clr, rain_fallspeed_cld,              &
-                          snow_fallspeed_clr, snow_fallspeed_cld
+                          snow_fallspeed_clr, snow_fallspeed_cld,              &
+                          ok_ice_sedim, ice_fallspeed_sedim
 
 
@@ -1364,4 +1368,5 @@
 REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
 REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qice_sedim     !--quantity of sedimentated ice [kg/kg]
 
 REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
@@ -1388,4 +1393,5 @@
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
 
 
@@ -1428,4 +1434,7 @@
 REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
 REAL :: Eff_rain_ice
+
+!--Ice sedimentation
+REAL :: dz, qice_sedim_tmp
 
 
@@ -1444,4 +1453,5 @@
 dqrfreez(:) = 0.
 dqsfreez(:) = 0.
+IF ( ok_ice_sedim ) dqised(:) = 0.
 
 
@@ -1559,5 +1569,5 @@
         !--tau for snow depends on the ice fraction in mixed-phase clouds     
         tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
 
         expo_auto_rain = cld_expo_con
@@ -1570,5 +1580,5 @@
         !--tau for snow depends on the ice fraction in mixed-phase clouds     
         tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
 
         expo_auto_rain = cld_expo_lsc
@@ -1903,4 +1913,24 @@
 
 
+  !---------------------------------------------------------
+  !--                  ICE SEDIMENTATION
+  !---------------------------------------------------------
+  IF ( ok_ice_sedim ) THEN
+    qice_sedim_tmp = 0.
+
+    IF ( temp(i) .LE. temp_nowater ) THEN
+      rho = pplay(i) / temp(i) / RD
+      dz = ( paprsdn(i) - paprsup(i) ) / RG / rho
+      qice_sedim_tmp = qice(i) * MIN(1., ice_fallspeed_sedim * dtime / dz)
+      !--Add tendencies
+      qice(i) = qice(i) - qice_sedim_tmp
+    ENDIF
+
+    !--Diagnostic tendencies
+    dqised(i) = ( qice_sedim(i) - qice_sedim_tmp ) / dtime
+    !--Save quantity that will be sedimentated in the layer below
+    qice_sedim(i) = qice_sedim_tmp
+  ENDIF
+
 
   !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
@@ -1915,5 +1945,4 @@
   !--Note that this is physically different than what is proposed in LTP thesis.
   precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
-  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
 
   !--Calculate outputs
Index: LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90	(revision 5589)
@@ -683,4 +683,14 @@
       REAL, SAVE, ALLOCATABLE :: dcf_avi(:,:), dqi_avi(:,:), dqvc_avi(:,:)
       !$OMP THREADPRIVATE(dcf_avi, dqi_avi, dqvc_avi)
+      REAL, SAVE, ALLOCATABLE :: cldfra_nocont(:,:), cldtau_nocont(:,:), cldemi_nocont(:,:)
+      !$OMP THREADPRIVATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
+      REAL, SAVE, ALLOCATABLE :: cldh_nocont(:), contcov(:)
+      !$OMP THREADPRIVATE(cldh_nocont, contcov)
+      REAL, SAVE, ALLOCATABLE :: fiwp_nocont(:), fiwc_nocont(:,:), ref_ice_nocont(:,:)
+      !$OMP THREADPRIVATE(fiwp_nocont, fiwc_nocont, ref_ice_nocont)
+      REAL, SAVE, ALLOCATABLE :: topsw_nocont(:), solsw_nocont(:)
+      !$OMP THREADPRIVATE(topsw_nocont, solsw_nocont)
+      REAL, SAVE, ALLOCATABLE :: toplw_nocont(:), sollw_nocont(:)
+      !$OMP THREADPRIVATE(toplw_nocont, sollw_nocont)
 
 !-- LSCP - mixed phase clouds variables
@@ -717,10 +727,10 @@
       REAL, SAVE, ALLOCATABLE :: dqsfreez(:,:)
       !$OMP THREADPRIVATE(dqsfreez)
-      REAL, SAVE, ALLOCATABLE :: dcfres(:,:)
-      !$OMP THREADPRIVATE(dcfres)
-      REAL, SAVE, ALLOCATABLE :: dqsres(:,:)
-      !$OMP THREADPRIVATE(dqsres)
-      REAL, SAVE, ALLOCATABLE :: dqvcres(:,:)
-      !$OMP THREADPRIVATE(dqvcres)
+      REAL, SAVE, ALLOCATABLE :: dqised(:,:)
+      !$OMP THREADPRIVATE(dqised)
+      REAL, SAVE, ALLOCATABLE :: dcfsed(:,:)
+      !$OMP THREADPRIVATE(dcfsed)
+      REAL, SAVE, ALLOCATABLE :: dqvcsed(:,:)
+      !$OMP THREADPRIVATE(dqvcsed)
 
 ! variables for stratospheric aerosol
@@ -1238,4 +1248,9 @@
       ALLOCATE(contfra(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))
+      ALLOCATE(cldh_nocont(klon), contcov(klon))
+      ALLOCATE(fiwp_nocont(klon), fiwc_nocont(klon,klev), ref_ice_nocont(klon,klev))
+      ALLOCATE(topsw_nocont(klon), solsw_nocont(klon))
+      ALLOCATE(toplw_nocont(klon), sollw_nocont(klon))
 
 !-- LSCP - POPRECIP variables
@@ -1244,5 +1259,5 @@
       ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
       ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
-      ALLOCATE(dcfres(klon,klev), dqsres(klon,klev), dqvcres(klon,klev))
+      ALLOCATE(dqised(klon,klev), dcfsed(klon,klev), dqvcsed(klon,klev))
 
 IF (CPPKEY_STRATAER) THEN
@@ -1646,4 +1661,7 @@
       DEALLOCATE(contfra, dcontfra_cir)
       DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
+      DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
+      DEALLOCATE(cldh_nocont, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)
+      DEALLOCATE(topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
 
 !-- LSCP - POPRECIP variables
@@ -1652,5 +1670,5 @@
       DEALLOCATE(dqrauto, dqrcol, dqrmelt, dqrfreez)
       DEALLOCATE(dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
-      DEALLOCATE(dcfres, dqsres, dqvcres)
+      DEALLOCATE(dqised, dcfsed, dqvcsed)
 
 IF (CPPKEY_STRATAER) THEN
Index: LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5589)
@@ -1611,10 +1611,10 @@
   TYPE(ctrl_out), SAVE :: o_dqsfreez = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'dqsfreez', 'LS snow tendency due to freezing', 'kg/kg/s', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_dcfres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
-    'dcfres', 'LS cloud fraction tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_dqsres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
-    'dqsres', 'LS snow tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_dqvcres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
-    'dqvcres', 'LS cloud water vapor tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqised = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqised', 'Ice water content tendency due to sedimentation', 'kg/kg/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dcfsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dcfsed', 'Cloud fraction tendency due to sedimentation', 's-1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqvcsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqvcsed', 'Cloud water vapor tendency due to sedimentation', 'kg/kg/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))
@@ -2214,4 +2214,30 @@
   TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'flighth2o', 'Aviation H2O flight emission', 'kg H2O/s/m^3', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_cldfra_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'cldfra_nocont', 'Cloud fraction w/o contrails', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_cldtau_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'cldtau_nocont', 'Cloud optical thickness w/o contrails', '1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_cldemi_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'cldemi_nocont', 'Cloud optical emissivity w/o contrails', '1', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_cldh_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'cldh_nocont', 'High-level cloudiness w/o contrails', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_contcov = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'contcov', 'Total contrails cover', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_iwp_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'iwp_nocont', 'Cloud ice water path w/o contrails', 'kg/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_iwc_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'iwc_nocont', 'Cloud ice water content seen by radiation w/o contrails', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_ref_ice_nocont = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'ref_ice_nocont', 'Effective radius of ice crystals w/o contrails', 'microns', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_tops_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'tops_nocont', 'Solar rad. at TOA w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_topl_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'topl_nocont', 'IR rad. at TOA w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_sols_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'sols_nocont', 'Solar rad. at surf. w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_soll_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'soll_nocont', 'IR rad. at surface w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_nettop_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'nettop_nocont', 'Net dn radiatif flux at TOA w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
 
 !!!!!!!!!!!!! Sorties niveaux standards de pression NMC 
Index: LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90	(revision 5589)
@@ -145,5 +145,5 @@
          o_qrainlsc, o_qsnowlsc, o_dqreva, o_dqrauto, o_dqrcol, o_dqrmelt, o_dqrfreez, &
          o_dqssub, o_dqsauto, o_dqsagg, o_dqsrim, o_dqsmelt, o_dqsfreez, &
-         o_dcfres, o_dqsres, o_dqvcres, &
+         o_dqised, o_dcfsed, o_dqvcsed, &
          o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, &
          o_dqsphy, o_dqsphy2d, o_dqbsphy, o_dqbsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
@@ -232,4 +232,7 @@
          o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
          o_dcontfracir, o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
+         o_cldfra_nocont, o_cldtau_nocont, o_cldemi_nocont, o_cldh_nocont, &
+         o_contcov, o_iwp_nocont, o_iwc_nocont, o_ref_ice_nocont, &
+         o_tops_nocont, o_topl_nocont, o_sols_nocont, o_soll_nocont, &
 !--interactive CO2
          o_flx_co2_ocean, o_flx_co2_ocean_cor, &
@@ -269,5 +272,5 @@
          o_SAD_sulfate, o_reff_sulfate, o_sulfmmr, o_nd_mode, o_sulfmmr_mode
 
-    USE lmdz_lscp_ini, ONLY: ok_poprecip, ok_precip_resuspension
+    USE lmdz_lscp_ini, ONLY: ok_poprecip, ok_ice_sedim
 
     USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
@@ -358,4 +361,7 @@
          Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
          dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
+         cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
+         contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
+         topsw_nocont, toplw_nocont, solsw_nocont, sollw_nocont, &
          alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
          alp_bl_stat, alp_bl_fluct_tke, slab_wfbils, &
@@ -384,5 +390,5 @@
          dqrauto,dqrcol,dqrmelt,dqrfreez, &
          dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
-         dcfres, dqsres, dqvcres, &
+         dqised, dcfsed, dqvcsed, &
          d_t_dyn,  & 
          d_q_dyn,  d_ql_dyn, d_qs_dyn, d_qbs_dyn,  &
@@ -2095,8 +2101,8 @@
            CALL histwrite_phy(o_dqsfreez, dqsfreez)
            CALL histwrite_phy(o_dqsrim, dqsrim)
-           IF ( ok_precip_resuspension ) THEN
-            CALL histwrite_phy(o_dcfres, dcfres)
-            CALL histwrite_phy(o_dqsres, dqsres)
-            CALL histwrite_phy(o_dqvcres, dqvcres)
+           IF ( ok_ice_sedim ) THEN
+            CALL histwrite_phy(o_dqised, dqised)
+            CALL histwrite_phy(o_dcfsed, dcfsed)
+            CALL histwrite_phy(o_dqvcsed, dqvcsed)
            ENDIF
            ELSE
@@ -2154,4 +2160,22 @@
          CALL histwrite_phy(o_dqiavi, dqi_avi)
          CALL histwrite_phy(o_dqvcavi, dqvc_avi)
+         CALL histwrite_phy(o_cldfra_nocont, cldfra_nocont)
+         CALL histwrite_phy(o_cldtau_nocont, cldtau_nocont)
+         CALL histwrite_phy(o_cldemi_nocont, cldemi_nocont)
+         CALL histwrite_phy(o_cldh_nocont, cldh_nocont)
+         CALL histwrite_phy(o_contcov, contcov)
+         CALL histwrite_phy(o_iwp_nocont, fiwp_nocont)
+         CALL histwrite_phy(o_iwc_nocont, fiwc_nocont)
+         CALL histwrite_phy(o_ref_ice_nocont, ref_ice_nocont)
+         IF (ok_rad_contrail) THEN
+           IF (vars_defined) zx_tmp_fi2d = topsw_nocont * swradcorr
+           CALL histwrite_phy(o_tops_nocont, topsw_nocont)
+           CALL histwrite_phy(o_topl_nocont, toplw_nocont)
+           IF (vars_defined) zx_tmp_fi2d = topsw_nocont * swradcorr - toplw_nocont
+           CALL histwrite_phy(o_nettop_nocont, zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d = solsw_nocont * swradcorr
+           CALL histwrite_phy(o_sols_nocont, solsw_nocont)
+           CALL histwrite_phy(o_soll_nocont, sollw_nocont)
+         ENDIF
        ENDIF
        IF (ok_plane_h2o) THEN
Index: LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90	(revision 5589)
@@ -187,5 +187,5 @@
        ! proprecip
        qraindiag, qsnowdiag, &
-       dqreva, dqssub, dcfres, dqsres, dqvcres, &
+       dqreva, dqssub, dqised, dcfsed, dqvcsed, &
        dqrauto,dqrcol,dqrmelt,dqrfreez, &
        dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
@@ -331,4 +331,7 @@
        contfra, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
        dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, &
+       cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
+       contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
+       topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont, &
        !
        stratomask, &
@@ -1226,4 +1229,6 @@
     ! "CRF off"
     REAL, dimension(klon, klev) :: cldfrarad   ! fraction nuageuse
+    !--AB contrails
+    REAL, dimension(klon, klev) :: cldfrarad_nocont   ! fraction nuageuse without contrails
 
     REAL :: calday, zxsnow_dummy(klon)
@@ -1430,4 +1435,10 @@
           WRITE (lunout, *) ' ok_plane_contrail=y requires 6 H2O tracers ', &
               '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c, H2O_a) but nqo=', nqo, '. Might as well stop here.'
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+
+       IF (ok_rad_contrail.AND..NOT.ok_plane_contrail) THEN
+          WRITE (lunout, *) ' ok_rad_contrail=y requires ok_plane_contrail=y '
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -2111,4 +2122,7 @@
        !=============================================================
 
+       !--Read the aviation emissions
+       IF ( ok_plane_h2o .OR. ok_plane_contrail ) CALL read_aviation_emissions(klon, klev)
+
        IF (using_xios) THEN
          ! Get "missing_val" value from XML files (from temperature variable)
@@ -2123,7 +2137,4 @@
          IF (is_omp_master) CALL xios_get_field_attr("temp",default_value=missing_val)
          CALL bcast_omp(missing_val)
-
-         !--Read the aviation emissions
-         IF ( ok_plane_h2o .OR. ok_plane_contrail ) CALL read_aviation_emissions(klon, klev)
        !
        ! Now we activate some double radiation call flags only if some
@@ -3924,5 +3935,5 @@
          qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
          dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez, &
-         dcfres, dqsres, dqvcres)
+         dqised, dcfsed, dqvcsed)
 
 
@@ -4470,5 +4481,8 @@
                ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
                reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
-               zfice, dNovrN, ptconv, rnebcon, clwcon, rcont_seri)
+               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)
 
        !
@@ -4479,4 +4493,6 @@
        cldemirad   = cldemi
        cldfrarad   = cldfra
+       !--AB contrails
+       IF (ok_rad_contrail) cldfrarad_nocont = cldfra_nocont
 
        !
@@ -4505,4 +4521,13 @@
              ENDDO
           ENDDO
+
+          !--AB contrails
+          IF ( ok_rad_contrail ) THEN
+            DO k=1, klev
+              DO i=1, klon
+                cldfrarad_nocont(i,k) = cldfra_nocont(i,k) * beta(i,k)
+              ENDDO
+            ENDDO
+          ENDIF
           !
        ELSE
@@ -4533,4 +4558,13 @@
              ENDDO
           ENDDO
+
+          !--AB contrails
+          IF ( ok_rad_contrail ) THEN
+            DO k=1, klev
+              DO i=1, klon
+                cldfrarad_nocont(i,k) = cldfra_nocont(i,k) * beta(i,k)
+              ENDDO
+            ENDDO
+          ENDIF
        !
        ENDIF
@@ -4649,5 +4683,8 @@
                ZLWFT0_i, ZFLDN0, ZFLUP0, &
                ZSWFT0_i, ZFSDN0, ZFSUP0, &
-               cloud_cover_sw)
+               cloud_cover_sw, &
+               !--AB contrails radiative effects
+               cldfrarad_nocont, fiwc_nocont, ref_ice_nocont, &
+               topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
 
           !lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
@@ -4660,4 +4697,9 @@
                         sollwdown(:))
           cool = cool + betalwoff * (cool0 - cool)
+          !--AB contrails
+          IF ( ok_rad_contrail ) THEN
+            toplw_nocont = toplw_nocont + betalwoff * (toplw0 - toplw_nocont)
+            sollw_nocont = sollw_nocont + betalwoff * (sollw0 - sollw_nocont)
+          ENDIF
 
           IF (.NOT. using_xios) THEN
@@ -4727,5 +4769,8 @@
                      ZLWFT0_i, ZFLDN0, ZFLUP0, &
                      ZSWFT0_i, ZFSDN0, ZFSUP0, &
-                     cloud_cover_sw)
+                     cloud_cover_sw, &
+                     !--AB contrails radiative effects
+                     cldfrarad_nocont, fiwc_nocont, ref_ice_nocont, &
+                     topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
           ENDIF !ok_4xCO2atm
 
Index: LMDZ6/branches/contrails/libf/phylmd/radlwsw_m.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/radlwsw_m.F90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/radlwsw_m.F90	(revision 5589)
@@ -45,5 +45,8 @@
        ZLWFT0_i, ZFLDN0, ZFLUP0, &
        ZSWFT0_i, ZFSDN0, ZFSUP0, &
-       cloud_cover_sw)
+       cloud_cover_sw, &
+       !--AB contrails radiative effects
+       cldfra_nocont, fiwc_nocont, ref_ice_nocont, &
+       topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
 
     ! Modules necessaires
@@ -279,4 +282,12 @@
     REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
     REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
+    !--AB contrails radiative effects
+    REAL, DIMENSION(klon,klev), INTENT(in)  :: cldfra_nocont
+    REAL, DIMENSION(klon,klev), INTENT(in)  :: fiwc_nocont
+    REAL, DIMENSION(klon,klev), INTENT(in)  :: ref_ice_nocont
+    REAL, DIMENSION(klon),      INTENT(out) :: topsw_nocont
+    REAL, DIMENSION(klon),      INTENT(out) :: solsw_nocont
+    REAL, DIMENSION(klon),      INTENT(out) :: toplw_nocont
+    REAL, DIMENSION(klon),      INTENT(out) :: sollw_nocont
 
     ! Local variables
@@ -467,4 +478,12 @@
     REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
     REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
+    !--AB contrails radiative effects
+    REAL(KIND=8) cldfra_nocont_i(klon,klev)
+    REAL(KIND=8) fiwc_nocont_i(klon,klev)
+    REAL(KIND=8) ref_ice_nocont_i(klon,klev)
+    REAL(KIND=8) ZTOPSWNOCONT(klon)
+    REAL(KIND=8) ZSOLSWNOCONT(klon)
+    REAL(KIND=8) ZTOPLWNOCONT(klon)
+    REAL(KIND=8) ZSOLLWNOCONT(klon)
     ! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
     !      REAL(KIND=8) RSUN(3,2)
@@ -900,4 +919,10 @@
              ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
              ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
+             IF (ok_rad_contrail) THEN
+               !--AB contrails radiative effects
+               cldfra_nocont_i(1:klon,k)  = cldfra_nocont(1:klon,klev+1-k)
+               fiwc_nocont_i(1:klon,k)    = fiwc_nocont(1:klon,klev+1-k)
+               ref_ice_nocont_i(1:klon,k) = ref_ice_nocont(1:klon,klev+1-k)
+             ENDIF
           ENDDO
           DO k=1,kflev
@@ -980,5 +1005,8 @@
                ZLWADAERO, & !--NL
                volmip_solsw, flag_volc_surfstrat, & !--VOLMIP
-               ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
+               ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback, & ! flags aerosols
+               !--AB contrails radiative effect
+               ok_rad_contrail, cldfra_nocont_i, fiwc_nocont_i, ref_ice_nocont_i, &
+               ZTOPSWNOCONT, ZSOLSWNOCONT, ZTOPLWNOCONT, ZSOLLWNOCONT)
 
           !--OB diagnostics
@@ -1080,4 +1108,10 @@
           ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
           ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
+
+          IF (ok_rad_contrail) THEN
+            !--AB contrails radiative effect
+            ZTOPSWNOCONT(:) = ZTOPSWNOCONT(:) * fract(:)
+            ZSOLSWNOCONT(:) = ZSOLSWNOCONT(:) * fract(:)
+          ENDIF
 
           ! ---------
@@ -1655,4 +1689,13 @@
           ENDDO
        ENDIF
+       !--AB radiative effect from contrails
+       IF (ok_rad_contrail) THEN
+         DO i = 1, kdlon
+           topsw_nocont(iof+i) = ZTOPSWNOCONT(i)
+           solsw_nocont(iof+i) = ZSOLSWNOCONT(i)
+           toplw_nocont(iof+i) = ZTOPLWNOCONT(i)
+           sollw_nocont(iof+i) = ZSOLLWNOCONT(i)
+         ENDDO
+       ENDIF
        DO k = 1, kflev
           DO i = 1, kdlon
Index: LMDZ6/branches/contrails/libf/phylmd/rrtm/recmwf_aero.F90
===================================================================
--- LMDZ6/branches/contrails/libf/phylmd/rrtm/recmwf_aero.F90	(revision 5585)
+++ LMDZ6/branches/contrails/libf/phylmd/rrtm/recmwf_aero.F90	(revision 5589)
@@ -40,5 +40,8 @@
      !..end
      & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
-     & flag_aer_feedback)
+     & flag_aer_feedback, &
+     !--AB contrails radiative effect
+     & ok_rad_contrail, PCLFR_NOCONT, PQIWP_NOCONT, PREF_ICE_NOCONT, &
+     & PTOPSWNOCONT, PSOLSWNOCONT, PTOPLWNOCONT, PSOLLWNOCONT)
   !--fin
 
@@ -266,4 +269,11 @@
   REAL(KIND=JPRB)   ,INTENT(OUT)   :: volmip_solsw(KPROMA) ! SW clear sky in the case of VOLMIP
   INTEGER, INTENT(IN)              :: flag_volc_surfstrat !--VOlMIP Modif
+  !--AB contrails radiative effect
+  LOGICAL           ,INTENT(IN)    :: ok_rad_contrail
+  REAL(KIND=JPRB)   ,INTENT(IN)    :: PCLFR_NOCONT(KPROMA,KLEV) 
+  REAL(KIND=JPRB)   ,INTENT(IN)    :: PQIWP_NOCONT(KPROMA,KLEV)
+  REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_ICE_NOCONT(KPROMA,KLEV)
+  REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWNOCONT(KPROMA), PSOLSWNOCONT(KPROMA) ! No contrails experiment forcing at TOA and surface (SW)
+  REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPLWNOCONT(KPROMA), PSOLLWNOCONT(KPROMA) ! No contrails experiment forcing at TOA and surface (LW)
 
   !     ==== COMPUTED IN RADITE ===
@@ -346,4 +356,11 @@
   REAL(KIND=JPRB) ::  LWUP0_AERO(KPROMA,KLEV+1,5)
   REAL(KIND=JPRB) ::  LWDN0_AERO(KPROMA,KLEV+1,5)
+  !--AB contrails radiative effect
+  REAL(KIND=JPRB) :: ZRCLC_NOCONT(KPROMA,KLEV), ZQIWP_NOCONT(KPROMA,KLEV)
+  REAL(KIND=JPRB) :: PREF_LIQ_NOCONT(KPROMA,KLEV)
+  REAL(KIND=JPRB) :: PPIZA_NOCONT(KPROMA,KLEV,NSW)
+  REAL(KIND=JPRB) :: PCGA_NOCONT(KPROMA,KLEV,NSW)
+  REAL(KIND=JPRB) :: PTAU_NOCONT(KPROMA,KLEV,NSW)
+  REAL(KIND=JPRB) :: PTAU_LW_NOCONT(KPROMA,KLEV,NLW)
 
 #include "radlsw.intfb.h"
@@ -673,4 +690,77 @@
   ENDIF ! .not. AEROSOLFEEDBACK_ACTIVE
 
+
+  !--Double call to radiative routine for contrails
+  !--The calculation are done again WITHOUT contrails
+  IF (ok_rad_contrail) THEN
+
+     !--The same base case is used
+     !--NB. here we could use pointers...
+     IF ( flag_aerosol .EQ. 0 ) THEN
+        PREF_LIQ_NOCONT(:,:) = PREF_LIQ_PI(:,:)
+        PPIZA_NOCONT(:,:,:) = PPIZA_ZERO(:,:,:)
+        PCGA_NOCONT(:,:,:) = PCGA_ZERO(:,:,:)
+        PTAU_NOCONT(:,:,:) = PTAU_ZERO(:,:,:)
+        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_ZERO(:,:,:)
+     ELSEIF ( .not. ok_ade .AND. .not. ok_aie ) THEN
+        PREF_LIQ_NOCONT(:,:) = PREF_LIQ_PI(:,:)
+        PPIZA_NOCONT(:,:,:) = PPIZA_NAT(:,:,:)
+        PCGA_NOCONT(:,:,:) = PCGA_NAT(:,:,:)
+        PTAU_NOCONT(:,:,:) = PTAU_NAT(:,:,:)
+        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_NAT(:,:,:)
+     ELSEIF ( .not. ok_ade .AND. ok_aie ) THEN
+        PREF_LIQ_NOCONT(:,:) = PREF_LIQ(:,:)
+        PPIZA_NOCONT(:,:,:) = PPIZA_NAT(:,:,:)
+        PCGA_NOCONT(:,:,:) = PCGA_NAT(:,:,:)
+        PTAU_NOCONT(:,:,:) = PTAU_NAT(:,:,:)
+        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_NAT(:,:,:)
+     ELSEIF ( ok_ade .AND. .not. ok_aie ) THEN
+        PREF_LIQ_NOCONT(:,:) = PREF_LIQ_PI(:,:)
+        PPIZA_NOCONT(:,:,:) = PPIZA_TOT(:,:,:)
+        PCGA_NOCONT(:,:,:) = PCGA_TOT(:,:,:)
+        PTAU_NOCONT(:,:,:) = PTAU_TOT(:,:,:)
+        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_TOT(:,:,:)
+     ELSEIF ( ok_ade .AND. ok_aie ) THEN
+        PREF_LIQ_NOCONT(:,:) = PREF_LIQ(:,:)
+        PPIZA_NOCONT(:,:,:) = PPIZA_TOT(:,:,:)
+        PCGA_NOCONT(:,:,:) = PCGA_TOT(:,:,:)
+        PTAU_NOCONT(:,:,:) = PTAU_TOT(:,:,:)
+        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_TOT(:,:,:)
+     ENDIF
+
+     DO JK=1,KLEV
+         DO JL=IBEG,IEND
+             ZRCLC_NOCONT(JL,JK)=MAX( 0.0_JPRB ,MIN( 1.0_JPRB ,PCLFR_NOCONT(JL,JK)))
+             IF (ZRCLC_NOCONT(JL,JK) > REPCLC) THEN
+                 ZQIWP_NOCONT(JL,JK)=PQIWP_NOCONT(JL,JK)
+             ELSE
+                 ZQIWP_NOCONT(JL,JK)=REPH2O*ZRCLC_NOCONT(JL,JK)
+             ENDIF
+         ENDDO
+     ENDDO
+
+     CALL RADLSW (&
+          & IBEG  , IEND   , KPROMA  , KLEV  , KMODE , NAER,&
+          & ZRII0 ,&
+          & ZRAER , PALBD  , PALBP   , PAPRS , ZRPR  ,&
+          & ZCCNL , ZCCNO  ,&
+          & PCCO2 , ZRCLC_NOCONT  , PDP     , PEMIS , ZEMIW ,PSLM    , ZRMU0 , ZPQO3,&
+          & ZQ    , ZQIWP_NOCONT  , ZQLWP   , ZQS   , ZQRAIN,ZQRAINT ,&
+          & PTH   , ZRTI   , PTS     , ZNBAS , ZNTOP ,&
+          & PREF_LIQ_NOCONT, PREF_ICE_NOCONT,&
+          & ZEMIT , ZFCT   , ZFLT    , ZFCS    , ZFLS  ,&
+          & ZFRSOD, ZSUDU  , ZUVDF   , ZPARF   , ZPARCF, ZTINCF, PSFSWDIR,&
+          & PSFSWDIF,PFSDNN, PFSDNV  ,&  
+          & LRDUST,PPIZA_NOCONT,PCGA_NOCONT,PTAU_NOCONT,PTAU_LW_NOCONT,PFLUX,PFLUC,&
+          & PFSDN , PFSUP  , PFSCDN  , PFSCUP )
+
+     ! save budgets LW and SW at TOA and surface
+     PTOPSWNOCONT(:) = PFSDN(:,KLEV+1) - PFSUP(:,KLEV+1)
+     PSOLSWNOCONT(:) = PFSDN(:,1) - PFSUP(:,1)
+     PTOPLWNOCONT(:) = PFLUX(:,1,KLEV+1) + PFLUX(:,2,KLEV+1)
+     PSOLLWNOCONT(:) = - PFLUX(:,1,1) - PFLUX(:,2,1)
+
+  ENDIF
+
   !*         4.2     TRANSFORM FLUXES TO MODEL HISTORICAL VARIABLES
 
