Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3474)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3479)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 MODULE etat0phys
 !
@@ -111,5 +114,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
@@ -134,5 +137,6 @@
                    iflag_cldcon,                                        &
                    iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
-                   ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple,     &
+                   ok_ade, ok_aie, ok_volcan, ok_alw, ok_cdnc,          &
+                   aerosol_couple,                                      &
                    chemistry_couple, flag_aerosol, flag_aerosol_strat,  & 
                    flag_aer_feedback,                                   &
Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 3479)
@@ -17,5 +17,5 @@
        iflag_cld_th, &
        iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
-       ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple, &
+       ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, chemistry_couple, &
        flag_aerosol, flag_aerosol_strat, flag_aer_feedback, new_aod, &
        flag_bc_internal_mixture, bl95_b0, bl95_b1,&
@@ -63,4 +63,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)
     !
 
@@ -71,5 +72,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
@@ -93,5 +94,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
@@ -2310,4 +2321,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
@@ -2656,4 +2668,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/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3479)
@@ -1380,4 +1380,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/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3479)
@@ -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)
@@ -216,4 +216,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, &
@@ -392,5 +395,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
@@ -1445,4 +1448,15 @@
        ENDIF
 #endif
+       !NL
+       IF (ok_volcan .AND. ok_ade) THEN
+          DO k=1, klev
+             IF (vars_defined) zx_tmp_fi3d(:,k)=heat_volc(:,k)*swradcorr(:)
+          ENDDO
+          CALL histwrite_phy(o_heat_volc, zx_tmp_fi3d)
+          DO k=1, klev
+             IF (vars_defined) zx_tmp_fi3d(:,k)=cool_volc(:,k)*swradcorr(:)
+          ENDDO
+          CALL histwrite_phy(o_cool_volc, zx_tmp_fi3d)
+       ENDIF
        IF (ok_ade) THEN
           IF (vars_defined) zx_tmp_fi2d(:)=topswad_aero*swradcorr
Index: LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3479)
@@ -314,4 +314,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)
@@ -324,4 +326,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)
@@ -575,4 +581,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))
@@ -718,4 +725,5 @@
       deallocate(heat, heat0) 
       deallocate(cool, cool0)
+      deallocate(heat_volc, cool_volc) 
       deallocate(topsw, toplw)
       deallocate(sollwdown, sollwdownclr)
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3479)
@@ -331,4 +331,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.)
@@ -1233,5 +1234,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, &
+            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
+            chemistry_couple, &
             flag_aerosol, flag_aerosol_strat, flag_aer_feedback, new_aod, &
             flag_bc_internal_mixture, bl95_b0, bl95_b1, &
@@ -1612,5 +1614,5 @@
        CALL 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,           &
@@ -3984,5 +3986,6 @@
                t_seri,q_seri,wo, &
                cldfrarad, cldemirad, cldtaurad, &
-               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
+               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, &
+               flag_aerosol, &
                flag_aerosol_strat, flag_aer_feedback, &
                tau_aero, piz_aero, cg_aero, &
@@ -3995,4 +3998,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, &
@@ -4069,5 +4073,6 @@
                      t_seri,q_seri,wo, &
                      cldfrarad, cldemirad, cldtaurad, &
-                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
+                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, &
+                     flag_aerosol, &
                      flag_aerosol_strat, flag_aer_feedback, &
                      tau_aero, piz_aero, cg_aero, &
@@ -4080,4 +4085,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, &
@@ -4940,5 +4946,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/trunk/libf/phylmd/radlwsw_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/radlwsw_m.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/radlwsw_m.F90	(revision 3479)
@@ -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, flag_aer_feedback, &
    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)
@@ -120,4 +122,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:
@@ -193,4 +198,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)
@@ -228,4 +234,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)
@@ -294,4 +301,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)
@@ -308,4 +316,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
@@ -430,4 +441,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.
@@ -573,8 +586,9 @@
       ENDDO
       DO k = 1, kflev
-      DO i = 1, kdlon
-      zcool(i,k)=0.
-      zcool0(i,k)=0.
-      ENDDO
+         DO i = 1, kdlon
+            zcool(i,k)=0.
+            zcool_volc(i,k)=0. !NL
+            zcool0(i,k)=0.
+         ENDDO
       ENDDO
       DO i = 1, kdlon
@@ -599,27 +613,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
@@ -869,8 +885,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, flag_aer_feedback) ! flags aerosols
+         ZLWADAERO, & !--NL
+         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
             
 !        print *,'RADLWSW: apres RECMWF'
@@ -951,4 +969,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
@@ -1031,4 +1053,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)
@@ -1140,4 +1166,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/trunk/libf/phylmd/regr_pr_time_av_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90	(revision 3479)
@@ -251,4 +251,5 @@
   END IF
   !$OMP END MASTER
+  !$OMP BARRIER
   IF(lfirst) THEN
     lfirst=.FALSE.;       CALL bcast(lfirst)
Index: LMDZ6/trunk/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3479)
@@ -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/trunk/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3479)
@@ -358,16 +358,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/trunk/libf/phylmd/rrtm/recmwf_aero.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/rrtm/recmwf_aero.F90	(revision 3474)
+++ LMDZ6/trunk/libf/phylmd/rrtm/recmwf_aero.F90	(revision 3479)
@@ -30,10 +30,13 @@
  & 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,flag_aer_feedback)
+ & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
+ & flag_aer_feedback)
 !--fin
 
@@ -82,4 +85,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)
@@ -213,4 +217,5 @@
 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
@@ -221,4 +226,7 @@
 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
@@ -811,4 +819,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
@@ -831,4 +842,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
@@ -845,4 +859,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
@@ -865,5 +882,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
@@ -879,5 +899,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))
@@ -899,5 +922,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))
@@ -913,5 +939,8 @@
      PSOLSWAD0AERO(:) = 0.0
      PTOPSWAD0AERO(:) = 0.0
-
+     IF(ok_volcan) THEN
+        PSWADAERO(:,:)  = 0.0 !--NL
+     ENDIF
+     
 ! indirect anthropogenic forcing
      PSOLSWAIAERO(:) = 0.0 
@@ -933,5 +962,8 @@
      PSOLLWAD0AERO(:) = 0.0
      PTOPLWAD0AERO(:) = 0.0
-
+     IF(ok_volcan) THEN
+        PLWADAERO(:,:)  = 0.0 !--NL
+     ENDIF
+     
 ! LW indirect anthropogenic forcing
      PSOLLWAIAERO(:) = 0.0 
