Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3425)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/radlwsw_m.F90	(revision 3426)
@@ -722,6 +722,6 @@
 !
 !--OB
-!--aerosol TOT  - anthropogenic+natural
-!--aerosol NAT  - natural only
+!--aerosol TOT  - anthropogenic+natural - index 2
+!--aerosol NAT  - natural only          - index 1
 !
       DO i = 1, kdlon
@@ -743,6 +743,6 @@
 !
 !--C. Kleinschmitt
-!--aerosol TOT  - anthropogenic+natural
-!--aerosol NAT  - natural only
+!--aerosol TOT  - anthropogenic+natural - index 2
+!--aerosol NAT  - natural only          - index 1
 !
       DO i = 1, kdlon
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 3425)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3426)
@@ -795,75 +795,57 @@
          cg_allaer(i,k,2,inu)=MIN(MAX(cg_allaer(i,k,2,inu),0.0),1.0)
 
-!--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)
+         IF (.NOT. ok_volcan) THEN
+!
+!--this is the default case
+!--in this case, index 1 of tau_allaer contains natural aerosols only
+!--because the objective is to perform the double radiation call with and without anthropogenic aerosols 
+!
+           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
+           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)
+           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_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)
+!
+!--this is the case for VOLMIP
+!--in this case index 1 of tau_allaer contains all (natural+anthropogenic) aerosols (same as index 2 above)
+!--but stratospheric aerosols will not be added in rrtm/readaerosolstrato2 as
+!--the objective is to have the double radiation call with and without stratospheric aerosols
+!
+           tau_allaer(i,k,1,inu)=tau_allaer(i,k,2,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
+           piz_allaer(i,k,1,inu)=piz_allaer(i,k,2,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)
+           cg_allaer(i,k,1,inu) =cg_allaer(i,k,2,inu)
+!
          ENDIF
         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 3425)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3426)
@@ -291,9 +291,12 @@
     ENDDO
 
-    IF ( .not. ok_volcan) THEN
+    IF (.NOT. ok_volcan) THEN
+!
+!--this is the default case 
+!--stratospheric aerosols are added to both index 2 and 1 for double radiation calls
 !--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 
+!--strat aerosols are added to index 2 : natural and anthropogenic aerosols for 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) ) /              &
@@ -304,7 +307,5 @@
                                     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
-! BEWARE: if you change stings here, you'll also need to change them a few 
-! lines down after the ELSE
+!--strat aerosols are added to index 1 : natural aerosols only for 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) ) /              &
@@ -317,10 +318,13 @@
     ENDWHERE
     ENDDO
-!--no stratospheric aerosol in index 1 for these tests
+!
     ELSE
+!
+!--this is the VOLMIP case
+!--stratospheric aerosols are only added to index 2 in this case 
 !--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 
+!--strat aerosols are added to index 2 : natural and anthropogenic aerosols for 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) ) /              &
@@ -331,8 +335,4 @@
                                     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
@@ -352,5 +352,7 @@
     ENDDO
 
-    IF ( .not. ok_volcan) THEN
+    IF (.NOT. ok_volcan) THEN
+!--this is the default case 
+!--stratospheric aerosols are added to both index 2 and 1
     DO band=1, NLW
       WHERE (stratomask.GT.0.999999)
@@ -359,9 +361,14 @@
       ENDWHERE
     ENDDO
+!
     ELSE
+!
+!--this is the VOLMIP case
     DO band=1, NLW
       WHERE (stratomask.GT.0.999999)
+!--stratospheric aerosols are only added to index 2
         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
+!--stratospheric aerosols are not added to index 1 
+!--and we copy index 2 in index 1 because we want the same dust aerosol LW properties as above
         tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,2,band)
       ENDWHERE
