Index: /LMDZ6/branches/IPSLCM6.0.15/DefLists/CMIP6_ping_atmos.xml
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/DefLists/CMIP6_ping_atmos.xml	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/DefLists/CMIP6_ping_atmos.xml	(revision 3408)
@@ -155,5 +155,5 @@
    <field id="CMIP6_loadss"        field_ref="loadss"           /> <!-- P1 (kg m-2) atmosphere_mass_content_of_seasalt_dry_aerosol : Load of Seasalt -->
    <field id="CMIP6_longitude"     field_ref="dummy_not_provided"         /> <!-- P1 (degrees_east) longitude : Longitude -->
-   <field id="CMIP6_lwsffluxaero"  field_ref="dummy_not_provided"         /> <!-- P2 (W m-2) longwave__flux__due_to_volcanic_aerosols_at_the_surface : downwelling longwave  flux  due to volcanic aerosols at the surface to be diagnosed through double radiation call -->
+   <field id="CMIP6_lwsffluxaero"  field_ref="toplwad"         /> <!-- P2 (W m-2) longwave__flux__due_to_volcanic_aerosols_at_the_surface : downwelling longwave  flux  due to volcanic aerosols at the surface to be diagnosed through double radiation call -->
    <field id="CMIP6_lwsrfasdust"   field_ref="dummy_not_provided"         /> <!-- P1 (W m-2) surface_instantaneous_longwave_forcing_due_to_dust : All-sky Surface Longwave radiative flux due to Dust -->
    <field id="CMIP6_lwsrfcsdust"   field_ref="dummy_not_provided"         /> <!-- P1 (W m-2) surface_instantaneous_longwave_forcing_due_to_dust_in_clearsky : Clear-sky Surface Longwave radiative flux due to Dust -->
@@ -161,5 +161,5 @@
    <field id="CMIP6_lwtoacsaer"    field_ref="toplwad0"         /> <!-- P1 (W m-2) toa_instantaneous_longwave_forcing : Clear-Sky LW-RF Aerosols at TOA -->
    <field id="CMIP6_lwtoacsdust"   field_ref="dummy_not_provided"         /> <!-- P1 (W m-2) toa_instantaneous_longwave_forcing_due_to_dust_in_clearsky : Clear-sky TOA Longwave radiative flux due to Dust -->
-   <field id="CMIP6_lwtoafluxaerocs" field_ref="dummy_not_provided"       /> <!-- P1 (W m-2) longwave_flux_due_to_volcanic_aerosols_at_TOA_under_clear_sky : downwelling longwave flux due to volcanic aerosols at TOA under clear sky to be diagnosed through double radiation call -->
+   <field id="CMIP6_lwtoafluxaerocs" field_ref="toplwad0"       /> <!-- P1 (W m-2) longwave_flux_due_to_volcanic_aerosols_at_TOA_under_clear_sky : downwelling longwave flux due to volcanic aerosols at TOA under clear sky to be diagnosed through double radiation call -->
    <field id="CMIP6_mc"            field_ref="mc"               /> <!-- P1 (kg m-2 s-1) atmosphere_net_upward_convective_mass_flux : The net mass flux should represent the difference between the updraft and downdraft components.  The flux is computed as the mass divided by the area of the grid cell. -->
    <field id="CMIP6_mcd"           field_ref="dnwd"              > (dnwd-dnwd0) &gt; 0 ? (dnwd-dnwd0) : 0 </field> <!-- P2 (kg m-2 s-1) atmosphere_downdraft_convective_mass_flux : Calculated as the convective mass flux divided by the area of the whole grid cell (not just the area of the cloud). -->
@@ -281,10 +281,10 @@
    <field id="CMIP6_snwc"          field_ref="dummy_not_provided"         /> <!-- P1 (kg m-2) canopy_snow_amount : canopy_snow_amount -->
    <field id="CMIP6_solbnd"        field_ref="solbnd"         />  <!-- P1 (W m-2) solar_irradiance : Top-of-Atmosphere Solar Insolation for each band -->
-   <field id="CMIP6_swsffluxaero"  field_ref="dummy_not_provided"         /> <!-- P2 (W m-2) shortwave__flux_due_to_volcanic_aerosols_at__the_surface : downwelling shortwave  flux due to volcanic aerosols at  the surface to be diagnosed through double radiation call -->
+   <field id="CMIP6_swsffluxaero"  field_ref="topswad"         /> <!-- P2 (W m-2) shortwave__flux_due_to_volcanic_aerosols_at__the_surface : downwelling shortwave  flux due to volcanic aerosols at  the surface to be diagnosed through double radiation call -->
    <field id="CMIP6_swsrfasdust"   field_ref="dummy_not_provided"         /> <!-- P1 (W m-2) tendency_of_all_sky_surface_shortwave_flux_due_to_dust_ambient_aerosol_particles : All-sky Surface Shortwave radiative flux due to Dust -->
    <field id="CMIP6_swsrfcsdust"   field_ref="dummy_not_provided"         /> <!-- P1 (W m-2) tendency_of_clear_sky_surface_shortwave_flux_due_to_dust_ambient_aerosol_particles : Clear-sky Surface Shortwave radiative flux due to Dust -->
    <field id="CMIP6_swtoaasdust"   field_ref="dummy_not_provided"         /> <!-- P1 (W m-2) toa_instantaneous_shortwave_forcing : all sky sw-rf dust at toa -->
    <field id="CMIP6_swtoacsdust"   field_ref="dummy_not_provided"         /> <!-- P1 (W m-2) toa_instantaneous_shortwave_forcing : clear sky sw-rf dust at toa -->
-   <field id="CMIP6_swtoafluxaerocs" field_ref="dummy_not_provided"       /> <!-- P1 (W m-2) shortwave_flux_due_to_volcanic_aerosols_at_TOA_under_clear_sky : downwelling shortwave flux due to volcanic aerosols at TOA under clear sky to be diagnosed through double radiation call -->
+   <field id="CMIP6_swtoafluxaerocs" field_ref="topswad0"       /> <!-- P1 (W m-2) shortwave_flux_due_to_volcanic_aerosols_at_TOA_under_clear_sky : downwelling shortwave flux due to volcanic aerosols at TOA under clear sky to be diagnosed through double radiation call -->
    <field id="CMIP6_sza"           field_ref="sza"              /> <!-- P1 (degree) solar_zenith_angle : solar zenith angle -->
    <field id="CMIP6_t2"            field_ref="temp"> temp*temp   </field> <!-- P2 (K2) square_of_air_temperature : square_of_air_temperature -->
@@ -363,6 +363,6 @@
    <field id="CMIP6_zhalf"         field_ref="zhalf"            /> <!-- P2 (m) height_above_reference_ellipsoid : This is actual height above mean sea level, not geopotential height.  This is actual height above mean sea level, not geopotential height.  Includes both the top of the model atmosphere and surface levels. -->
    <field id="CMIP6_zmla"          field_ref="s_pblh"           /> <!-- P1 (m) atmosphere_boundary_layer_thickness : Height of Boundary Layer -->
-   <field id="CMIP6_zmlwaero"      field_ref="dummy_not_provided"        /> <!-- P1 (K s-1) longwave_heating_rate_due_to_volcanic_aerosols : longwave heating rate due to volcanic aerosols to be diagnosed through double radiation call, zonal average values required -->
-   <field id="CMIP6_zmswaero"      field_ref="dummy_not_provided"        /> <!-- P1 (K s-1) shortwave_heating_rate_due_to_volcanic_aerosols : shortwave heating rate due to volcanic aerosols to be diagnosed through double radiation call, zonal average values required -->
+   <field id="CMIP6_zmlwaero"      field_ref="cool_volc"        /> <!-- P1 (K s-1) longwave_heating_rate_due_to_volcanic_aerosols : longwave heating rate due to volcanic aerosols to be diagnosed through double radiation call, zonal average values required -->
+   <field id="CMIP6_zmswaero"      field_ref="heat_volc"        /> <!-- P1 (K s-1) shortwave_heating_rate_due_to_volcanic_aerosols : shortwave heating rate due to volcanic aerosols to be diagnosed through double radiation call, zonal average values required -->
    <field id="CMIP6_zmtnt"         field_ref="dtphy"            /> <!-- P1 (K s-1) tendency_of_air_temperature_due_to_diabatic_processes : Zonal Mean Diabatic Heating Rates -->
    <field id="CMIP6_ap" field_ref="Ahyb" /><!-- Ap hybrid coordinate array for level interfaces -->
Index: /LMDZ6/branches/IPSLCM6.0.15/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/DefLists/field_def_lmdz.xml	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/DefLists/field_def_lmdz.xml	(revision 3408)
@@ -702,4 +702,6 @@
         <field id="rldcs4co2"    long_name="Downwelling CS LW 4xCO2 atmosphere"    unit="W/m2" />
         <field id="stratomask"    long_name="Stratospheric fraction"    unit="-" />
+        <field id="heat_volc"    long_name="SW heating rate from volcano"    unit="K/s" />
+        <field id="cool_volc"    long_name="LW cooling rate from volcano"    unit="K/s" />
         <field id="pvap"    long_name="pvap intermediary variable" unit="-">pres*ovap*461.5 / (287.04*(1.+ (10.9491/18.0153)*ovap)) </field>
     </field_group>
Index: /LMDZ6/branches/IPSLCM6.0.15/DefLists/file_def_histday_lmdz.xml
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/DefLists/file_def_histday_lmdz.xml	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/DefLists/file_def_histday_lmdz.xml	(revision 3408)
@@ -337,12 +337,12 @@
                 <field field_ref="z0m" level="10" />
                 <field field_ref="z0h" level="10" />
-                <field field_ref="topswad" level="10" />
-                <field field_ref="topswad0" level="10" />
+                <field field_ref="topswad" level="1" />
+                <field field_ref="topswad0" level="1" />
                 <field field_ref="topswai" level="10" />
                 <field field_ref="solswad" level="10" />
                 <field field_ref="solswad0" level="10" />
                 <field field_ref="solswai" level="10" />
-                <field field_ref="toplwad" level="10" />
-                <field field_ref="toplwad0" level="10" />
+                <field field_ref="toplwad" level="1" />
+                <field field_ref="toplwad0" level="1" />
                 <field field_ref="toplwai" level="10" />
                 <field field_ref="sollwad" level="10" />
@@ -628,4 +628,7 @@
                 <field field_ref="rsdcs4co2" level="10" />
                 <field field_ref="rldcs4co2" level="10" />
+                <field field_ref="stratomask" level="1" />
+                <field field_ref="cool_volc" level="1" />
+                <field field_ref="heat_volc" level="1" />
             </field_group>
         </file>
Index: /LMDZ6/branches/IPSLCM6.0.15/DefLists/file_def_histmth_lmdz.xml
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/DefLists/file_def_histmth_lmdz.xml	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/DefLists/file_def_histmth_lmdz.xml	(revision 3408)
@@ -372,12 +372,12 @@
                 <field field_ref="z0m" level="10" />
                 <field field_ref="z0h" level="10" />
-                <field field_ref="topswad" level="10" />
-                <field field_ref="topswad0" level="10" />
+                <field field_ref="topswad" level="1" />
+                <field field_ref="topswad0" level="1" />
                 <field field_ref="topswai" level="10" />
                 <field field_ref="solswad" level="10" />
                 <field field_ref="solswad0" level="10" />
                 <field field_ref="solswai" level="10" />
-                <field field_ref="toplwad" level="10" />
-                <field field_ref="toplwad0" level="10" />
+                <field field_ref="toplwad" level="1" />
+                <field field_ref="toplwad0" level="1" />
                 <field field_ref="toplwai" level="10" />
                 <field field_ref="sollwad" level="10" />
@@ -677,4 +677,7 @@
                 <field field_ref="rsdcs4co2" level="10" />
                 <field field_ref="rldcs4co2" level="10" />
+                <field field_ref="stratomask" level="1" />
+                <field field_ref="cool_volc" level="1" />
+                <field field_ref="heat_volc" level="1" />
             </field_group>
 
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3408)
@@ -110,5 +110,5 @@
   INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
   REAL    :: ratqsbas, ratqshaut, tau_ratqs
-  LOGICAL :: ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
+  LOGICAL :: ok_ade, ok_aie, ok_volcan, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
   INTEGER :: flag_aerosol
   INTEGER :: flag_aerosol_strat
@@ -126,14 +126,14 @@
 ! Physics configuration
 !*******************************************************************************
-  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
-                   callstats,                                           &
-                   solarlong0,seuil_inversion,                          &
-                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
-                   iflag_cldcon,                                        &
-                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
-                   ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple,     &
-                   chemistry_couple, flag_aerosol, flag_aerosol_strat,  &
-                   new_aod, flag_bc_internal_mixture, bl95_b0, bl95_b1, &
-                   read_climoz, alp_offset)
+  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,       &
+                   callstats,                                             &
+                   solarlong0,seuil_inversion,                            &
+                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,        &
+                   iflag_cldcon,                                          &
+                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,              &
+                   ok_ade, ok_aie, ok_volcan, ok_alw, ok_cdnc,            &
+                   aerosol_couple, chemistry_couple, flag_aerosol,        &
+                   flag_aerosol_strat, new_aod, flag_bc_internal_mixture, &
+                   bl95_b0, bl95_b1, read_climoz, alp_offset)
   CALL phys_state_var_init(read_climoz)
 
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90	(revision 3408)
@@ -17,6 +17,6 @@
        iflag_cld_th, &
        iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
-       ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple, &
-       flag_aerosol, flag_aerosol_strat, new_aod, &
+       ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
+       chemistry_couple, flag_aerosol, flag_aerosol_strat, new_aod, &
        flag_bc_internal_mixture, bl95_b0, bl95_b1,&
        read_climoz, &
@@ -66,4 +66,5 @@
     ! flag_bc_internal_mixture : use BC internal mixture if true
     ! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc 
+    ! ok_volcan: activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
     !
 
@@ -75,5 +76,5 @@
     LOGICAL              :: ok_LES
     LOGICAL              :: callstats
-    LOGICAL              :: ok_ade, ok_aie, ok_alw, ok_cdnc
+    LOGICAL              :: ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan
     LOGICAL              :: aerosol_couple, chemistry_couple
     INTEGER              :: flag_aerosol
@@ -96,5 +97,5 @@
     LOGICAL, SAVE       :: ok_LES_omp   
     LOGICAL, SAVE       :: callstats_omp
-    LOGICAL, SAVE       :: ok_ade_omp, ok_aie_omp, ok_alw_omp, ok_cdnc_omp
+    LOGICAL, SAVE       :: ok_ade_omp, ok_aie_omp, ok_alw_omp, ok_cdnc_omp, ok_volcan_omp
     LOGICAL, SAVE       :: aerosol_couple_omp, chemistry_couple_omp
     INTEGER, SAVE       :: flag_aerosol_omp
@@ -396,4 +397,14 @@
     ok_cdnc_omp = .FALSE.
     CALL getin('ok_cdnc', ok_cdnc_omp)
+
+    !
+    !Config Key  = ok_volcan
+    !Config Desc = ok to generate volcanic diags
+    !Config Def  = .FALSE.
+    !Config Help = Used in radlwsw_m.F
+    !
+    ok_volcan_omp = .FALSE.
+    CALL getin('ok_volcan', ok_volcan_omp)
+
     !
     !Config Key  = aerosol_couple
@@ -2282,4 +2293,5 @@
     ok_alw = ok_alw_omp
     ok_cdnc = ok_cdnc_omp
+    ok_volcan = ok_volcan_omp
     aerosol_couple = aerosol_couple_omp
     chemistry_couple = chemistry_couple_omp
@@ -2613,4 +2625,5 @@
     write(lunout,*)' pmagic = ',pmagic
     write(lunout,*)' ok_ade = ',ok_ade
+    write(lunout,*)' ok_volcan = ',ok_volcan
     write(lunout,*)' ok_aie = ',ok_aie
     write(lunout,*)' ok_alw = ',ok_alw
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3408)
@@ -1375,4 +1375,8 @@
   TYPE(ctrl_out), SAVE :: o_temp = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
     'temp', 'Air temperature', 'K', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_heat_volc = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'heat_volc', 'SW heating rate due to volcano', 'K/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_cool_volc = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'cool_volc', 'LW cooling rate due to volcano', 'K/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_theta = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
     'theta', 'Potential air temperature', 'K', (/ ('', i=1, 10) /))
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_output_write_mod.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_output_write_mod.F90	(revision 3408)
@@ -17,5 +17,5 @@
   SUBROUTINE phys_output_write(itap, pdtphys, paprs, pphis, &
        pplay, lmax_th, aerosol_couple,         &
-       ok_ade, ok_aie, ivap, iliq, isol, new_aod, ok_sync, &
+       ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, ok_sync, &
        ptconv, read_climoz, clevSTD, ptconvth, &
        d_u, d_t, qx, d_qx, zmasse, flag_aerosol, flag_aerosol_strat, ok_cdnc)
@@ -214,4 +214,7 @@
          o_vsed_aer, o_tau_strat_1020, o_ext_strat_1020, o_f_r_wet
 #endif
+    
+    USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
+    USE phys_state_var_mod, ONLY: heat_volc, cool_volc !NL
 
     USE phys_state_var_mod, ONLY: pctsrf, rain_fall, snow_fall, &
@@ -388,5 +391,5 @@
     INTEGER, DIMENSION(klon) :: lmax_th
     LOGICAL :: aerosol_couple, ok_sync
-    LOGICAL :: ok_ade, ok_aie, new_aod
+    LOGICAL :: ok_ade, ok_aie, ok_volcan, new_aod
     LOGICAL, DIMENSION(klon, klev) :: ptconv, ptconvth
     REAL :: pdtphys
@@ -1378,4 +1381,15 @@
        ENDIF
 #endif
+       !NL
+       IF (ok_volcan .AND. ok_ade) THEN
+          DO k=1, klev
+             zx_tmp_fi3d(:,k)=heat_volc(:,k)*swradcorr(:)
+          ENDDO
+          CALL histwrite_phy(o_heat_volc, zx_tmp_fi3d)
+          DO k=1, klev
+             zx_tmp_fi3d(:,k)=cool_volc(:,k)*swradcorr(:)
+          ENDDO
+          CALL histwrite_phy(o_cool_volc, zx_tmp_fi3d)
+       ENDIF
        IF (ok_ade) THEN
           CALL histwrite_phy(o_topswad, topswad_aero*swradcorr)
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_state_var_mod.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_state_var_mod.F90	(revision 3408)
@@ -311,4 +311,6 @@
 ! toplwdown : downward CS LW flux at TOA
 ! toplwdownclr : downward CS LW flux at TOA
+! heat_volc : chauffage solaire du au volcanisme
+! cool_volc : refroidissement infrarouge du au volcanisme
       REAL,ALLOCATABLE,SAVE :: clwcon0(:,:),rnebcon0(:,:)
 !$OMP THREADPRIVATE(clwcon0,rnebcon0)
@@ -321,4 +323,8 @@
       REAL,ALLOCATABLE,SAVE :: cool0(:,:)
 !$OMP THREADPRIVATE(cool0)
+      REAL,ALLOCATABLE,SAVE :: heat_volc(:,:)   
+!$OMP THREADPRIVATE(heat_volc)
+      REAL,ALLOCATABLE,SAVE :: cool_volc(:,:)
+!$OMP THREADPRIVATE(cool_volc)
       REAL,ALLOCATABLE,SAVE :: topsw(:), toplw(:)
 !$OMP THREADPRIVATE(topsw,toplw)
@@ -563,4 +569,5 @@
       ALLOCATE(heat(klon,klev), heat0(klon,klev)) 
       ALLOCATE(cool(klon,klev), cool0(klon,klev))
+      ALLOCATE(heat_volc(klon,klev), cool_volc(klon,klev)) 
       ALLOCATE(topsw(klon), toplw(klon))
       ALLOCATE(sollwdown(klon), sollwdownclr(klon))
@@ -698,4 +705,5 @@
       deallocate(heat, heat0) 
       deallocate(cool, cool0)
+      deallocate(heat_volc, cool_volc) 
       deallocate(topsw, toplw)
       deallocate(sollwdown, sollwdownclr)
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3408)
@@ -327,4 +327,5 @@
     include "dimpft.h"
     !======================================================================
+    LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques
     LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
     PARAMETER (ok_cvl=.TRUE.)
@@ -1219,6 +1220,6 @@
             fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
             iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
-            ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple,  chemistry_couple, &
-            flag_aerosol, flag_aerosol_strat, new_aod, &
+            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
+            chemistry_couple, flag_aerosol, flag_aerosol_strat, new_aod, &
             flag_bc_internal_mixture, bl95_b0, bl95_b1, &
                                 ! nv flags pour la convection et les
@@ -3879,6 +3880,6 @@
                t_seri,q_seri,wo, &
                cldfrarad, cldemirad, cldtaurad, &
-               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
-               flag_aerosol_strat, &
+               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
+               flag_aerosol, flag_aerosol_strat, &
                tau_aero, piz_aero, cg_aero, &
                tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & 
@@ -3890,4 +3891,5 @@
                ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
                heat,heat0,cool,cool0,albpla, &
+               heat_volc,cool_volc, &
                topsw,toplw,solsw,sollw, &
                sollwdown, &
@@ -3964,6 +3966,6 @@
                      t_seri,q_seri,wo, &
                      cldfrarad, cldemirad, cldtaurad, &
-                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
-                     flag_aerosol_strat, &
+                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
+                     flag_aerosol, flag_aerosol_strat, &
                      tau_aero, piz_aero, cg_aero, &
                      tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
@@ -3975,4 +3977,5 @@
                      ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
                      heatp,heat0p,coolp,cool0p,albplap, &
+                     heat_volc,cool_volc, &
                      topswp,toplwp,solswp,sollwp, &
                      sollwdownp, &
@@ -4823,5 +4826,5 @@
     CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
          pplay, lmax_th, aerosol_couple,                 &
-         ok_ade, ok_aie, ivap, iliq, isol, new_aod,      & 
+         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod,      & 
          ok_sync, ptconv, read_climoz, clevSTD,          &
          ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3408)
@@ -16,5 +16,5 @@
    t,q,wo,&
    cldfra, cldemi, cldtaupd,&
-   ok_ade, ok_aie, flag_aerosol,&
+   ok_ade, ok_aie, ok_volcan, flag_aerosol,&
    flag_aerosol_strat,&
    tau_aero, piz_aero, cg_aero,&
@@ -25,4 +25,5 @@
    ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
    heat,heat0,cool,cool0,albpla,&
+   heat_volc, cool_volc,&
    topsw,toplw,solsw,sollw,&
    sollwdown,&
@@ -100,4 +101,5 @@
   ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
   ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
+  ! ok_volcan-input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
   ! flag_aerosol-input-I- aerosol flag from 0 to 6
   ! flag_aerosol_strat-input-I- use stratospheric aerosols flag (0, 1, 2)
@@ -119,4 +121,7 @@
   ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
   ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
+  !
+  ! heat_volc-----output-R- echauffement atmospherique  du au forcage volcanique (visible) (K/s)
+  ! cool_volc-----output-R- refroidissement dans l'IR du au forcage volcanique (K/s)
   !
   ! ATTENTION: swai and swad have to be interpreted in the following manner:
@@ -192,4 +197,5 @@
 
   LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
+  LOGICAL, INTENT(in)  :: ok_volcan                                      ! produce volcanic diags (SW/LW heat flux and rate)
   LOGICAL              :: lldebug
   INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
@@ -226,4 +232,5 @@
   REAL,    INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
   REAL,    INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
+  REAL,    INTENT(out) :: heat_volc(KLON,KLEV), cool_volc(KLON,KLEV) !NL
   REAL,    INTENT(out) :: topsw(KLON), toplw(KLON)
   REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON)
@@ -292,4 +299,5 @@
   REAL(KIND=8) zheat(kdlon,kflev), zcool(kdlon,kflev)
   REAL(KIND=8) zheat0(kdlon,kflev), zcool0(kdlon,kflev)
+  REAL(KIND=8) zheat_volc(kdlon,kflev), zcool_volc(kdlon,kflev) !NL
   REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
   REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
@@ -306,4 +314,7 @@
   REAL(KIND=8) ztopswad0aero(kdlon), zsolswad0aero(kdlon)   ! Aerosol direct forcing at TOAand surface
   REAL(KIND=8) ztopswaiaero(kdlon), zsolswaiaero(kdlon)     ! dito, indirect
+!--NL
+  REAL(KIND=8) zswadaero(kdlon,kflev+1)                       ! SW Aerosol direct forcing
+  REAL(KIND=8) zlwadaero(kdlon,kflev+1)                       ! LW Aerosol direct forcing
 !-LW by CK
   REAL(KIND=8) ztoplwadaero(kdlon), zsollwadaero(kdlon)     ! LW Aerosol direct forcing at TOAand surface
@@ -413,4 +424,6 @@
       heat(i,k)=0.
       cool(i,k)=0.
+      heat_volc(i,k)=0. !NL
+      cool_volc(i,k)=0. !NL
       heat0(i,k)=0.
       cool0(i,k)=0.
@@ -544,27 +557,28 @@
 !--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------              
       DO k = 1, kflev+1
+         DO i = 1, kdlon
+!     print *,'RADLWSW: boucle mise a zero i k',i,k
+            ZFLUP(i,k)=0.
+            ZFLDN(i,k)=0.
+            ZFLUP0(i,k)=0.
+            ZFLDN0(i,k)=0.
+            ZLWFT0_i(i,k)=0.
+            ZFLUCUP_i(i,k)=0.
+            ZFLUCDWN_i(i,k)=0.
+         ENDDO
+      ENDDO
+      DO k = 1, kflev
+         DO i = 1, kdlon
+            zcool(i,k)=0.
+            zcool_volc(i,k)=0. !NL
+            zcool0(i,k)=0.
+         ENDDO
+      ENDDO
       DO i = 1, kdlon
-!     print *,'RADLWSW: boucle mise a zero i k',i,k
-      ZFLUP(i,k)=0.
-      ZFLDN(i,k)=0.
-      ZFLUP0(i,k)=0.
-      ZFLDN0(i,k)=0.
-      ZLWFT0_i(i,k)=0.
-      ZFLUCUP_i(i,k)=0.
-      ZFLUCDWN_i(i,k)=0.
-      ENDDO
-      ENDDO
-      DO k = 1, kflev
-      DO i = 1, kdlon
-      zcool(i,k)=0.
-      zcool0(i,k)=0.
-      ENDDO
-      ENDDO
-      DO i = 1, kdlon
-      ztoplw(i)=0.
-      zsollw(i)=0.
-      ztoplw0(i)=0.
-      zsollw0(i)=0.
-      zsollwdown(i)=0.
+         ztoplw(i)=0.
+         zsollw(i)=0.
+         ztoplw0(i)=0.
+         zsollw0(i)=0.
+         zsollwdown(i)=0.
       ENDDO
        ! Old radiation scheme, used for AR4 runs
@@ -582,27 +596,29 @@
 !----- Mise a zero des tableaux output du rayonnement SW-AR4 
       DO k = 1, kflev+1
-      DO i = 1, kdlon
-      ZFSUP(i,k)=0.
-      ZFSDN(i,k)=0.
-      ZFSUP0(i,k)=0.
-      ZFSDN0(i,k)=0.
-      ZFSUPC0(i,k)=0.
-      ZFSDNC0(i,k)=0.
-      ZFLUPC0(i,k)=0.
-      ZFLDNC0(i,k)=0.
-      ZSWFT0_i(i,k)=0.
-      ZFCUP_i(i,k)=0.
-      ZFCDWN_i(i,k)=0.
-      ZFCCUP_i(i,k)=0.
-      ZFCCDWN_i(i,k)=0.
-      ZFLCCUP_i(i,k)=0.
-      ZFLCCDWN_i(i,k)=0.
-      ENDDO
+         DO i = 1, kdlon
+            ZFSUP(i,k)=0.
+            ZFSDN(i,k)=0.
+            ZFSUP0(i,k)=0.
+            ZFSDN0(i,k)=0.
+            ZFSUPC0(i,k)=0.
+            ZFSDNC0(i,k)=0.
+            ZFLUPC0(i,k)=0.
+            ZFLDNC0(i,k)=0.
+            ZSWFT0_i(i,k)=0.
+            ZFCUP_i(i,k)=0.
+            ZFCDWN_i(i,k)=0.
+            ZFCCUP_i(i,k)=0.
+            ZFCCDWN_i(i,k)=0.
+            ZFLCCUP_i(i,k)=0.
+            ZFLCCDWN_i(i,k)=0.
+            zswadaero(i,k)=0. !--NL
+         ENDDO
       ENDDO
       DO k = 1, kflev
-      DO i = 1, kdlon
-      zheat(i,k)=0.
-      zheat0(i,k)=0.
-      ENDDO
+         DO i = 1, kdlon
+            zheat(i,k)=0.
+            zheat_volc(i,k)=0.
+            zheat0(i,k)=0.
+         ENDDO
       ENDDO
       DO i = 1, kdlon
@@ -852,8 +868,10 @@
          ZTOPSWAIAERO,ZSOLSWAIAERO, &
          ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
+         ZSWADAERO, & !--NL
          ZTOPLWADAERO,ZSOLLWADAERO,&  ! rajoute par C. Kleinscmitt pour LW diagnostics
          ZTOPLWAD0AERO,ZSOLLWAD0AERO,&
          ZTOPLWAIAERO,ZSOLLWAIAERO, &
-         ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat) ! flags aerosols
+         ZLWADAERO, & !--NL
+         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat) ! flags aerosols
             
 !        print *,'RADLWSW: apres RECMWF'
@@ -934,4 +952,8 @@
          ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
          ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
+         IF(ok_volcan) THEN
+            ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
+         ENDIF
+         
 !   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
 !   en sortie de radlsw.F90 - MPL 7.01.09
@@ -1014,4 +1036,8 @@
            zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
            zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
+           IF(ok_volcan) THEN
+              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RDAY*RG/RCPD/PDP(i,k) !NL
+              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RDAY*RG/RCPD/PDP(i,k) !NL
+           ENDIF
 !          print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
 !	   ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
@@ -1123,4 +1149,8 @@
         heat0(iof+i,k) = zheat0(i,k)/zznormcp
         cool0(iof+i,k) = zcool0(i,k)/zznormcp
+        IF(ok_volcan) THEN !NL
+           heat_volc(iof+i,k) = zheat_volc(i,k)/zznormcp
+           cool_volc(iof+i,k) = zcool_volc(i,k)/zznormcp
+        ENDIF
       ENDDO
     ENDDO
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3408)
@@ -794,40 +794,76 @@
          cg_allaer(i,k,2,inu)=MIN(MAX(cg_allaer(i,k,2,inu),0.0),1.0)
 
+! NL VOLC
 !--natural aerosol
-!--ASBCM aerosols take _pi value because of internal mixture option
-         tau_allaer(i,k,1,inu)=tau_ae_pi(i,k,id_ASSO4M_phy,inu)+tau_ae_pi(i,k,id_CSSO4M_phy,inu)+ &
-                               tau_ae_pi(i,k,id_ASBCM_phy,inu)+tau_ae_pi(i,k,id_AIBCM_phy,inu)+   &
-                               tau_ae_pi(i,k,id_ASPOMM_phy,inu)+tau_ae_pi(i,k,id_AIPOMM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_ASSSM_phy,inu)+tau_ae_pi(i,k,id_CSSSM_phy,inu)+   &
-                               tau_ae_pi(i,k,id_SSSSM_phy,inu)+ tau_ae_pi(i,k,id_CIDUSTM_phy,inu)
+! (same as upper but no volc aer in strat)
+         tau_allaer(i,k,1,inu)=tau_ae(i,k,id_ASSO4M_phy,inu)+tau_ae(i,k,id_CSSO4M_phy,inu)+ &
+                               tau_ae(i,k,id_ASBCM_phy,inu)+tau_ae(i,k,id_AIBCM_phy,inu)+   &
+                               tau_ae(i,k,id_ASPOMM_phy,inu)+tau_ae(i,k,id_AIPOMM_phy,inu)+ &
+                               tau_ae(i,k,id_ASSSM_phy,inu)+tau_ae(i,k,id_CSSSM_phy,inu)+   &
+                               tau_ae(i,k,id_SSSSM_phy,inu)+ tau_ae(i,k,id_CIDUSTM_phy,inu)
          tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),tau_min)
 
-         piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+   &
-                                tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+   &
-                                tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)+  &
-                                tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+     &
-                                tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+   &
-                                tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+   &
-                                tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+     &
-                                tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+     &
-                                tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+     &
-                                tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) &
+         piz_allaer(i,k,1,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+   &
+                                tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+   &
+                                tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+     &
+                                tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+     &
+                                tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+   &
+                                tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+   &
+                                tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+     &
+                                tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+     &
+                                tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+     &
+                                tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) &
                                 /tau_allaer(i,k,1,inu)
          piz_allaer(i,k,1,inu)=MIN(MAX(piz_allaer(i,k,1,inu),0.01),1.0)
          IF (tau_allaer(i,k,1,inu).LE.tau_min) piz_allaer(i,k,1,inu)=1.0
 
-         cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+    &
-                               tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+    &
-                               tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)*cg_ae_pi(i,k,id_ASBCM_phy,inu)+ &
-                               tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+       &
-                               tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+    &
-                               tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+    &
-                               tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+       &
-                               tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+       &
-                               tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+       &
-                               tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ &
+         cg_allaer(i,k,1,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ &
+                               tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ &
+                               tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+    &
+                               tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+    &
+                               tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ &
+                               tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ &
+                               tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+    &
+                               tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+    &
+                               tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+    &
+                               tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ &
                                (tau_allaer(i,k,1,inu)*piz_allaer(i,k,1,inu))
          cg_allaer(i,k,1,inu)=MIN(MAX(cg_allaer(i,k,1,inu),0.0),1.0)
 
+!--ASBCM aerosols take _pi value because of internal mixture option
+!         tau_allaer(i,k,1,inu)=tau_ae_pi(i,k,id_ASSO4M_phy,inu)+tau_ae_pi(i,k,id_CSSO4M_phy,inu)+ &
+!                               tau_ae_pi(i,k,id_ASBCM_phy,inu)+tau_ae_pi(i,k,id_AIBCM_phy,inu)+   &
+!                               tau_ae_pi(i,k,id_ASPOMM_phy,inu)+tau_ae_pi(i,k,id_AIPOMM_phy,inu)+ &
+!                               tau_ae_pi(i,k,id_ASSSM_phy,inu)+tau_ae_pi(i,k,id_CSSSM_phy,inu)+   &
+!                               tau_ae_pi(i,k,id_SSSSM_phy,inu)+ tau_ae_pi(i,k,id_CIDUSTM_phy,inu)
+!         tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),tau_min)
+!
+!         piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+   &
+!                                tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+   &
+!                                tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)+  &
+!                                tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+     &
+!                                tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+   &
+!                                tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+   &
+!                                tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+     &
+!                                tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+     &
+!                                tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+     &
+!                                tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) &
+!                                /tau_allaer(i,k,1,inu)
+!         piz_allaer(i,k,1,inu)=MIN(MAX(piz_allaer(i,k,1,inu),0.01),1.0)
+!         IF (tau_allaer(i,k,1,inu).LE.tau_min) piz_allaer(i,k,1,inu)=1.0
+!
+!         cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+    &
+!                               tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+    &
+!                               tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)*cg_ae_pi(i,k,id_ASBCM_phy,inu)+ &
+!                               tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+       &
+!                               tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+    &
+!                               tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+    &
+!                               tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+       &
+!                               tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+       &
+!                               tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+       &
+!                               tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ &
+!                               (tau_allaer(i,k,1,inu)*piz_allaer(i,k,1,inu))
+!         cg_allaer(i,k,1,inu)=MIN(MAX(cg_allaer(i,k,1,inu),0.0),1.0)
+! NL VOLC END
         ENDDO
       ENDDO
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3408)
@@ -303,16 +303,16 @@
         tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
 !--natural aerosols bands 1 to NSW
-        cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
-                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
-                                    MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
-                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
-        piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
-                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
-                                    MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
-        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
+!        cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
+!                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
+!                                    MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
+!                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
+!        piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
+!                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
+!                                    MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
+!        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
 !--no stratospheric aerosol in index 1 for these tests
-!        cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
-!        piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
-!        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
+        cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
+        piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
+        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
     ENDWHERE
     ENDDO
Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/recmwf_aero.F90
===================================================================
--- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/recmwf_aero.F90	(revision 3407)
+++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/recmwf_aero.F90	(revision 3408)
@@ -30,10 +30,12 @@
  & PTOPSWAIAERO,PSOLSWAIAERO,&
  & PTOPSWCFAERO,PSOLSWCFAERO,&
+ & PSWADAERO,& !--NL
 !--LW diagnostics CK
  & PTOPLWADAERO,PSOLLWADAERO,&
  & PTOPLWAD0AERO,PSOLLWAD0AERO,&
  & PTOPLWAIAERO,PSOLLWAIAERO,&
+ & PLWADAERO,& !--NL
 !..end
- & ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
+ & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat)
 !--fin
 
@@ -82,4 +84,5 @@
 ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
 ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
+! ok_volcan-input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
 ! flag_aerosol-input-I- aerosol flag from 0 to 7
 ! flag_aerosol_strat-input-I- use stratospheric aerosols flag (T/F)
@@ -212,11 +215,15 @@
 REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_ICE_PI(KPROMA,KLEV)
 LOGICAL, INTENT(in)  :: ok_ade, ok_aie         ! switches whether to use aerosol direct (indirect) effects or not
+LOGICAL, INTENT(in)  :: ok_volcan              ! produce volcanic diags (SW/LW heat flux and rate)
 INTEGER, INTENT(in)  :: flag_aerosol           ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
 LOGICAL, INTENT(in)  :: flag_aerosol_strat     ! use stratospheric aerosols
-REAL(KIND=JPRB)   ,INTENT(out)   :: PTOPSWADAERO(KPROMA), PSOLSWADAERO(KPROMA)       ! Aerosol direct forcing at TOA and surface
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWADAERO(KPROMA), PSOLSWADAERO(KPROMA)       ! Aerosol direct forcing at TOA and surface
 REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWAD0AERO(KPROMA), PSOLSWAD0AERO(KPROMA)     ! Aerosol direct forcing at TOA and surface
 REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWAIAERO(KPROMA), PSOLSWAIAERO(KPROMA)       ! ditto, indirect
 REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWCFAERO(KPROMA,3), PSOLSWCFAERO(KPROMA,3) !--do we keep this ?
 !--fin
+!--NL
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSWADAERO(KPROMA, KLEV+1)                        ! SW Aerosol direct forcing
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PLWADAERO(KPROMA, KLEV+1)                        ! LW Aerosol direct forcing
 !--CK
 REAL(KIND=JPRB)   ,INTENT(out)   :: PTOPLWADAERO(KPROMA), PSOLLWADAERO(KPROMA)       ! LW Aerosol direct forcing at TOA + surface
@@ -806,4 +813,7 @@
      PSOLSWAD0AERO(:) = (ZFSDN0_AERO(:,1,4)     -ZFSUP0_AERO(:,1,4))     -(ZFSDN0_AERO(:,1,2)     -ZFSUP0_AERO(:,1,2))
      PTOPSWAD0AERO(:) = (ZFSDN0_AERO(:,KLEV+1,4)-ZFSUP0_AERO(:,KLEV+1,4))-(ZFSDN0_AERO(:,KLEV+1,2)-ZFSUP0_AERO(:,KLEV+1,2))
+     IF(ok_volcan) THEN
+        PSWADAERO(:,:)  = (ZFSDN_AERO(:,:,4) -ZFSUP_AERO(:,:,4)) -(ZFSDN_AERO(:,:,2) -ZFSUP_AERO(:,:,2)) !--NL
+     ENDIF
 
 ! indirect anthropogenic forcing
@@ -826,4 +836,7 @@
      PSOLLWAD0AERO(:) = (-LWDN0_AERO(:,1,4)     -LWUP0_AERO(:,1,4))     -(-LWDN0_AERO(:,1,2)     -LWUP0_AERO(:,1,2))
      PTOPLWAD0AERO(:) = (-LWDN0_AERO(:,KLEV+1,4)-LWUP0_AERO(:,KLEV+1,4))-(-LWDN0_AERO(:,KLEV+1,2)-LWUP0_AERO(:,KLEV+1,2))
+     IF(ok_volcan) THEN
+        PLWADAERO(:,:)  = (-LWDN_AERO(:,:,4) -LWUP_AERO(:,:,4)) -(-LWDN_AERO(:,:,2) -LWUP_AERO(:,:,2)) !--NL
+     ENDIF
 
 ! LW indirect anthropogenic forcing
@@ -840,4 +853,7 @@
      PSOLSWAD0AERO(:) = (ZFSDN0_AERO(:,1,3)     -ZFSUP0_AERO(:,1,3))     -(ZFSDN0_AERO(:,1,1)     -ZFSUP0_AERO(:,1,1))
      PTOPSWAD0AERO(:) = (ZFSDN0_AERO(:,KLEV+1,3)-ZFSUP0_AERO(:,KLEV+1,3))-(ZFSDN0_AERO(:,KLEV+1,1)-ZFSUP0_AERO(:,KLEV+1,1))
+     IF(ok_volcan) THEN
+        PSWADAERO(:,:)  = (ZFSDN_AERO(:,:,3) -ZFSUP_AERO(:,:,3)) -(ZFSDN_AERO(:,:,1) -ZFSUP_AERO(:,:,1)) !--NL
+     ENDIF
 
 ! indirect anthropogenic forcing
@@ -860,5 +876,8 @@
      PSOLLWAD0AERO(:) = (-LWDN0_AERO(:,1,3)     -LWUP0_AERO(:,1,3))     -(-LWDN0_AERO(:,1,1)     -LWUP0_AERO(:,1,1))
      PTOPLWAD0AERO(:) = (-LWDN0_AERO(:,KLEV+1,3)-LWUP0_AERO(:,KLEV+1,3))-(-LWDN0_AERO(:,KLEV+1,1)-LWUP0_AERO(:,KLEV+1,1))
-
+     IF(ok_volcan) THEN
+        PLWADAERO(:,:)  = (-LWDN_AERO(:,:,3) -LWUP_AERO(:,:,3)) -(-LWDN_AERO(:,:,1) -LWUP_AERO(:,:,1)) !--NL
+     ENDIF
+     
 ! LW indirect anthropogenic forcing
      PSOLLWAIAERO(:) = 0.0
@@ -874,5 +893,8 @@
      PSOLSWAD0AERO(:) = 0.0
      PTOPSWAD0AERO(:) = 0.0
-
+     IF(ok_volcan) THEN
+        PSWADAERO(:,:)  = 0.0 !--NL
+     ENDIF
+     
 ! indirect anthropogenic forcing
      PSOLSWAIAERO(:) = (ZFSDN_AERO(:,1,2)     -ZFSUP_AERO(:,1,2))     -(ZFSDN_AERO(:,1,1)     -ZFSUP_AERO(:,1,1))
@@ -894,5 +916,8 @@
      PSOLLWAD0AERO(:) = 0.0
      PTOPLWAD0AERO(:) = 0.0
-
+     IF(ok_volcan) THEN
+        PLWADAERO(:,:)  = 0.0 !--NL
+     ENDIF
+     
 ! LW indirect anthropogenic forcing
      PSOLLWAIAERO(:) = (-LWDN_AERO(:,1,2)     -LWUP_AERO(:,1,2))     -(-LWDN_AERO(:,1,1)     -LWUP_AERO(:,1,1))
@@ -908,5 +933,8 @@
      PSOLSWAD0AERO(:) = 0.0
      PTOPSWAD0AERO(:) = 0.0
-
+     IF(ok_volcan) THEN
+        PSWADAERO(:,:)  = 0.0 !--NL
+     ENDIF
+     
 ! indirect anthropogenic forcing
      PSOLSWAIAERO(:) = 0.0 
@@ -928,5 +956,8 @@
      PSOLLWAD0AERO(:) = 0.0
      PTOPLWAD0AERO(:) = 0.0
-
+     IF(ok_volcan) THEN
+        PLWADAERO(:,:)  = 0.0 !--NL
+     ENDIF
+     
 ! LW indirect anthropogenic forcing
      PSOLLWAIAERO(:) = 0.0