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 3419)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_output_write_mod.F90	(revision 3425)
@@ -1388,5 +1388,5 @@
           CALL histwrite_phy(o_heat_volc, zx_tmp_fi3d)
           DO k=1, klev
-             zx_tmp_fi3d(:,k)=cool_volc(:,k)*swradcorr(:)
+             zx_tmp_fi3d(:,k)=cool_volc(:,k)
           ENDDO
           CALL histwrite_phy(o_cool_volc, zx_tmp_fi3d)
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3419)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3425)
@@ -1352,4 +1352,9 @@
 
        IF (type_trac == 'inca') THEN
+          ! jg : initialisation jusqu'au ces variables sont dans restart
+          ccm(:,:,:) = 0.
+          tau_aero(:,:,:,:) = 0.
+          piz_aero(:,:,:,:) = 0.
+          cg_aero(:,:,:,:) = 0.
 
           config_inca='none' ! default
@@ -3603,5 +3608,5 @@
 #else
                    !--climatologies or INCA aerosols
-                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, &
+                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
                         new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
                         pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
@@ -3693,5 +3698,5 @@
              CALL readaerosolstrato1_rrtm(debut)
             ELSEIF (flag_aerosol_strat.EQ.2) THEN
-             CALL readaerosolstrato2_rrtm(debut)
+             CALL readaerosolstrato2_rrtm(debut, ok_volcan)
             ELSE
              abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3419)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3425)
@@ -1037,6 +1037,6 @@
            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
+              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
+              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
            ENDIF
 !          print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
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 3419)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3425)
@@ -6,5 +6,5 @@
      tau_allaer, piz_allaer, &
      cg_allaer, m_allaer_pi, &
-     flag_aerosol, flag_bc_internal_mixture, zrho )
+     flag_aerosol, flag_bc_internal_mixture, zrho, ok_volcan )
 
   USE dimphy
@@ -32,4 +32,5 @@
   LOGICAL,                        INTENT(IN)  :: flag_bc_internal_mixture
   REAL, DIMENSION(klon,klev),     INTENT(IN)  :: zrho
+  LOGICAL,                        INTENT(IN)  :: ok_volcan ! volcanic diags
   !
   ! Output arguments:
@@ -794,15 +795,50 @@
          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
+         IF(.not. ok_volcan) THEN
+         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)
+         ELSE
 ! (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_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(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_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(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)+     &
@@ -815,8 +851,8 @@
                                 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(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,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(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)+    &
@@ -829,41 +865,6 @@
                                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
+            cg_allaer(i,k,1,inu)=MIN(MAX(cg_allaer(i,k,1,inu),0.0),1.0)
+         ENDIF
         ENDDO
       ENDDO
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90	(revision 3419)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90	(revision 3425)
@@ -1,5 +1,5 @@
 ! $Id$
 !
-SUBROUTINE readaerosol_optic_rrtm(debut, aerosol_couple, ok_alw,  &
+SUBROUTINE readaerosol_optic_rrtm(debut, aerosol_couple, ok_alw, ok_volcan, &
      new_aod, flag_aerosol, flag_bc_internal_mixture, itap, rjourvrai, &
      pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
@@ -32,4 +32,5 @@
   LOGICAL, INTENT(IN)                      :: aerosol_couple
   LOGICAL, INTENT(IN)                      :: ok_alw
+  LOGICAL, INTENT(IN)                      :: ok_volcan
   LOGICAL, INTENT(IN)                      :: new_aod
   INTEGER, INTENT(IN)                      :: flag_aerosol
@@ -313,5 +314,5 @@
        tau_aero, piz_aero, cg_aero, &
        m_allaer_pi, flag_aerosol,   &
-       flag_bc_internal_mixture, zrho ) 
+       flag_bc_internal_mixture, zrho, ok_volcan ) 
 
   ! aeropt_5wv only for validation and diagnostics
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3419)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3425)
@@ -2,5 +2,5 @@
 ! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
 !
-SUBROUTINE readaerosolstrato2_rrtm(debut)
+SUBROUTINE readaerosolstrato2_rrtm(debut, ok_volcan)
 
     USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, & 
@@ -29,4 +29,5 @@
 ! Variable input
     LOGICAL, INTENT(IN) ::  debut
+    LOGICAL, INTENT(IN) ::  ok_volcan !activate volcanic diags
 
 ! Variables locales
@@ -290,4 +291,5 @@
     ENDDO
 
+    IF ( .not. ok_volcan) THEN
 !--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
     DO band=1, NSW
@@ -303,18 +305,37 @@
         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)
+! BEWARE: if you change stings here, you'll also need to change them a few 
+! lines down after the ELSE
+        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)
+    ENDWHERE
+    ENDDO
 !--no stratospheric aerosol in index 1 for these tests
+    ELSE
+!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
+    DO band=1, NSW
+      WHERE (stratomask.GT.0.999999)
+!--anthropogenic aerosols bands 1 to NSW 
+        cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
+                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
+                                    MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
+                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
+        piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
+                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
+                                    MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
+        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)  = piz_aero_sw_rrtm(:,:,1,band)
         tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
-    ENDWHERE
-    ENDDO
+     ENDWHERE
+  ENDDO
+  ENDIF
 
 !--total vertical aod at 10 um
@@ -331,12 +352,20 @@
     ENDDO
 
+    IF ( .not. ok_volcan) THEN
     DO band=1, NLW
       WHERE (stratomask.GT.0.999999)
         tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
         tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
-!--no stratospheric aerosols in index 1 for these tests
-!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) 
       ENDWHERE
     ENDDO
+    ELSE
+    DO band=1, NLW
+      WHERE (stratomask.GT.0.999999)
+        tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
+        !--no stratospheric aerosols in index 1 for these tests
+        tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,2,band)
+      ENDWHERE
+    ENDDO
+    ENDIF
 
 !--default SSA value if there is no aerosol
