Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90	(revision 1237)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 SUBROUTINE AEROPT_2BANDS( &
      pdel, m_allaer, delt, RHcl, &
@@ -44,5 +47,5 @@
   REAL, DIMENSION(klon,klev,naero_tot,nbands) ::  cg_ae
   LOGICAL ::  soluble
-  INTEGER :: i, k, inu, m, mrfspecies
+  INTEGER :: i, k, ierr, inu, m, mrfspecies
   INTEGER :: spsol, spinsol, spss
   INTEGER :: RH_num
@@ -50,6 +53,6 @@
 
   INTEGER, PARAMETER :: nbre_RH=12
-  INTEGER, PARAMETER :: nbsol_compaer=3    ! 1- BC soluble; 2- POM soluble; 3- SO4.
-  INTEGER, PARAMETER :: nbinsol_compaer=3  ! 1- Dust; 2- BC insoluble; 3- POM insoluble
+  INTEGER, PARAMETER :: naero_soluble=7    ! 1- BC soluble; 2- POM soluble; 3- SO4.
+  INTEGER, PARAMETER :: naero_insoluble=3  ! 1- Dust; 2- BC insoluble; 3- POM insoluble
   LOGICAL, SAVE :: firstcall=.TRUE. 
 
@@ -70,5 +73,6 @@
 
 ! Coefficient optiques interpole sur le nombre de niveau du modele
-  REAL, DIMENSION(klev) :: A1_ASSSM_b1, A2_ASSSM_b1, A3_ASSSM_b1,&
+  REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
+          A1_ASSSM_b1, A2_ASSSM_b1, A3_ASSSM_b1,&
           B1_ASSSM_b1, B2_ASSSM_b1, C1_ASSSM_b1, C2_ASSSM_b1,&
           A1_CSSSM_b1, A2_CSSSM_b1, A3_CSSSM_b1,&
@@ -101,10 +105,10 @@
   ! Proprietes optiques
   !
-  REAL:: alpha_aers_2bands(nbre_RH,nbands,nbsol_compaer)   !--unit m2/g SO4
-  REAL:: alpha_aeri_2bands(nbands,nbinsol_compaer)
-  REAL:: cg_aers_2bands(nbre_RH,nbands,nbsol_compaer)      !--unit 
-  REAL:: cg_aeri_2bands(nbands,nbinsol_compaer)
-  REAL:: piz_aers_2bands(nbre_RH,nbands,nbsol_compaer)     !-- unit
-  REAL:: piz_aeri_2bands(nbands,nbinsol_compaer)           !-- unit
+  REAL:: alpha_aers_2bands(nbre_RH,nbands,naero_soluble)   !--unit m2/g SO4
+  REAL:: alpha_aeri_2bands(nbands,naero_insoluble)
+  REAL:: cg_aers_2bands(nbre_RH,nbands,naero_soluble)      !--unit 
+  REAL:: cg_aeri_2bands(nbands,naero_insoluble)
+  REAL:: piz_aers_2bands(nbre_RH,nbands,naero_soluble)     !-- unit
+  REAL:: piz_aeri_2bands(nbands,naero_insoluble)           !-- unit
 
 
@@ -354,6 +358,25 @@
        7.556,8.613,10.687,12.265,16.32,21.692, & 
        1.107,1.239,1.381,1.490,1.635,1.8030,   &
-       2.071,2.407,3.126,3.940,5.539,7.921/ 
-
+       2.071,2.407,3.126,3.940,5.539,7.921,    &
+                                ! sulfate coarse
+       4.681,5.062,5.460,5.798,6.224,6.733,    & 
+       7.556,8.613,10.687,12.265,16.32,21.692, & 
+       1.107,1.239,1.381,1.490,1.635,1.8030,   &
+       2.071,2.407,3.126,3.940,5.539,7.921,    &
+                                ! seasalt Super Coarse Soluble (SS)
+       0.5090,0.6554,0.7129,0.7767,0.8529,1.2728, &
+       1.3820,1.5792,1.9173,2.2002,2.7173,4.1487, &
+       0.5167,0.6613,0.7221,0.7868,0.8622,1.3027, &
+       1.4227,1.6317,1.9887,2.2883,2.8356,4.3453, &
+                                ! seasalt  Coarse Soluble (CS)
+       0.5090,0.6554,0.7129,0.7767,0.8529,1.2728, &
+       1.3820,1.5792,1.9173,2.2002,2.7173,4.1487, &
+       0.5167,0.6613,0.7221,0.7868,0.8622,1.3027, &
+       1.4227,1.6317,1.9887,2.2883,2.8356,4.3453, &
+                                ! seasalt  Accumulation Soluble (AS)
+       4.125, 4.674, 5.005, 5.434, 5.985, 10.006, &
+       11.175,13.376,17.264,20.540,26.604, 42.349,&
+       4.187, 3.939, 3.919, 3.937, 3.995,  5.078, &
+       5.511, 6.434, 8.317,10.152,14.024, 26.537/
 
   DATA alpha_aeri_2bands/  & 
@@ -364,5 +387,4 @@
        ! pom insoluble
        3.741,0.606/
-
 
   DATA cg_aers_2bands/ &
@@ -381,6 +403,25 @@
        .719, .733, .752, .760, .773, .786, &
        .544, .555, .565, .573, .583, .593, &
-       .610, .628, .655, .666, .692, .719/
-
+       .610, .628, .655, .666, .692, .719, &
+                                ! sulfate coarse
+       .658, .669, .680, .688, .698, .707, &
+       .719, .733, .752, .760, .773, .786, &
+       .544, .555, .565, .573, .583, .593, &
+       .610, .628, .655, .666, .692, .719, &
+                                ! seasalt Super Coarse soluble (SS)
+       .727, .747, .755, .761, .770, .788, &
+       .792, .799, .805, .809, .815, .826, &
+       .717, .738, .745, .752, .761, .779, &
+       .781, .786, .793, .797, .803, .813, &
+                                ! seasalt Coarse soluble (CS)
+       .727, .747, .755, .761, .770, .788, &
+       .792, .799, .805, .809, .815, .826, &
+       .717, .738, .745, .752, .761, .779, &
+       .781, .786, .793, .797, .803, .813, &
+                                ! Sesalt Accumulation Soluble (AS)
+       .727, .741, .748, .754, .761, .782, &
+       .787, .792, .797, .799, .801, .799, &
+       .606, .645, .658, .669, .681, .726, &
+       .734, .746, .761, .770, .782, .798/
 
   DATA cg_aeri_2bands/ &
@@ -407,6 +448,25 @@
        1.000,1.000,1.000,1.000,1.000,1.000, &
        .992, .988, .988, .987, .986, .985,  &
-       .985, .985, .984, .984, .984, .984 /
-
+       .985, .985, .984, .984, .984, .984,  &
+                                ! sulfate coarse
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       .992, .988, .988, .987, .986, .985,  &
+       .985, .985, .984, .984, .984, .984,  &
+                                ! seasalt Super Coarse Soluble (SS)
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       0.992,0.989,0.987,0.986,0.986,0.980, &
+       0.980,0.978,0.976,0.976,0.974,0.971, &
+                                ! seasalt Coarse soluble (CS)
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       0.992,0.989,0.987,0.986,0.986,0.980, &
+       0.980,0.978,0.976,0.976,0.974,0.971, &
+                                ! seasalt Accumulation Soluble (AS)
+       1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+       1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+       0.970, 0.975, 0.976, 0.977, 0.978, 0.982, &
+       0.982, 0.983, 0.984, 0.984, 0.985, 0.985/
 
   DATA piz_aeri_2bands/ &
@@ -418,8 +478,23 @@
        .966, .859/
 
-
 ! Interpolation des coefficients optiques de 19 niveaux vers le nombre des niveaux du model
   IF (firstcall) THEN
      firstcall=.FALSE.
+     
+     IF (.NOT. ALLOCATED(A1_ASSSM_b1)) THEN
+        ALLOCATE(A1_ASSSM_b1(klev),A2_ASSSM_b1(klev), A3_ASSSM_b1(klev),&
+          B1_ASSSM_b1(klev), B2_ASSSM_b1(klev), C1_ASSSM_b1(klev), C2_ASSSM_b1(klev),&
+          A1_CSSSM_b1(klev), A2_CSSSM_b1(klev), A3_CSSSM_b1(klev),&
+          B1_CSSSM_b1(klev), B2_CSSSM_b1(klev), C1_CSSSM_b1(klev), C2_CSSSM_b1(klev),&
+          A1_SSSSM_b1(klev), A2_SSSSM_b1(klev), A3_SSSSM_b1(klev),&
+          B1_SSSSM_b1(klev), B2_SSSSM_b1(klev), C1_SSSSM_b1(klev), C2_SSSSM_b1(klev),&
+          A1_ASSSM_b2(klev), A2_ASSSM_b2(klev), A3_ASSSM_b2(klev),&
+          B1_ASSSM_b2(klev), B2_ASSSM_b2(klev), C1_ASSSM_b2(klev), C2_ASSSM_b2(klev),&
+          A1_CSSSM_b2(klev), A2_CSSSM_b2(klev), A3_CSSSM_b2(klev),&
+          B1_CSSSM_b2(klev), B2_CSSSM_b2(klev), C1_CSSSM_b2(klev), C2_CSSSM_b2(klev),&
+          A1_SSSSM_b2(klev), A2_SSSSM_b2(klev), A3_SSSSM_b2(klev),&
+          B1_SSSSM_b2(klev), B2_SSSSM_b2(klev), C1_SSSSM_b2(klev), C2_SSSSM_b2(klev), stat=ierr)
+        IF (ierr /= 0) CALL abort_gcm('aeropt_2bands', 'pb in allocation 1',1)
+     END IF
      
 ! bande 1
@@ -478,16 +553,14 @@
   DO k=1, klev
      DO i=1, klon
-        IF (t_seri(i,k).EQ.0.) THEN
-           WRITE(lunout,*) 't_seri(i,k)=0 for i=',i,'k=',k
-           CALL abort_gcm('aeropt_2bands','t_seri=0',1)
-        END IF
-        IF (pplay(i,k).EQ.0.) THEN
-           WRITE(lunout,*) 'pplay(i,k)=0 for i=',i,'k=',k
-           CALL abort_gcm('aeropt_2bands','pplay=0',1)
-        END IF
-
-        zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
+!        IF (t_seri(i,k).EQ.0.) THEN
+!           WRITE(lunout,*) 't_seri(i,k)=0 for i=',i,'k=',k
+!           CALL abort_gcm('aeropt_2bands','t_seri=0',1)
+!        END IF
+!        IF (pplay(i,k).EQ.0.) THEN
+!           WRITE(lunout,*) 'pplay(i,k)=0 for i=',i,'k=',k
+!           CALL abort_gcm('aeropt_2bands','pplay=0',1)
+!        END IF
+        zrho=pplay(i,k)/t_seri(i,k)/RD                    ! kg/m3
         mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9
-
      ENDDO
   ENDDO
@@ -563,12 +636,17 @@
          spss=0
      ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN 
-         soluble=.TRUE.
-         spsol=2 
-         spss=0
-     ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR.  (aerosol_name(m).EQ.id_CSSO4M)) THEN
-         soluble=.TRUE.
-         spsol=3
-         spss=0
-         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+        soluble=.TRUE.
+        spsol=2 
+        spss=0
+     ELSEIF (aerosol_name(m).EQ.id_ASSO4M) THEN 
+        soluble=.TRUE.
+        spsol=3
+        spss=0
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+     ELSEIF  (aerosol_name(m).EQ.id_CSSO4M) THEN
+        soluble=.TRUE.
+        spsol=4
+        spss=0
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
      ELSEIF (aerosol_name(m).EQ.id_SSSSM) THEN 
          soluble=.TRUE.
@@ -604,17 +682,19 @@
      cg_ae2b_int(:,:,:)=0.
 
-     DO k=1, KLEV
-        DO i=1, KLON
-
-           rh=MIN(RHcl(i,k)*100.,RH_MAX)
-           RH_num = INT( rh/10. + 1.)
-
-           IF (rh.GT.85.) RH_num=10
-           IF (rh.GT.90.) RH_num=11
-           DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
-
-           DO inu=1,nbands
+     DO inu=1,nbands
+        DO k=1, KLEV
+           DO i=1, KLON
+
+              rh=MIN(RHcl(i,k)*100.,RH_MAX)
+              RH_num = INT( rh/10. + 1.)
+
+              IF (rh.GT.85.) RH_num=10
+              IF (rh.GT.90.) RH_num=11
+              DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
+
+!           DO inu=1,nbands
               IF (soluble) THEN
 
+              ! First optical parameters are computed for seasalt
                   IF (spss.NE.0) THEN
                       H=rh/100
@@ -664,11 +744,11 @@
                          DELTA* (cg_aers_2bands(RH_num+1,inu,spsol) - & 
                          cg_aers_2bands(RH_num,inu,spsol))
-                      
                   ENDIF
 
-                  tau_ae(i,k,aerosol_name(m),inu) = &
-                     mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*tau_ae2b_int(i,k,inu)*fac
-
-              ELSE
+                 tau_ae(i,k,aerosol_name(m),inu) = &
+                      mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt* &
+                                         tau_ae2b_int(i,k,inu)*fac
+
+              ELSE                                                    ! For all aerosol insoluble components
                  tau_ae2b_int(i,k,inu) = alpha_aeri_2bands(inu,spinsol)
                  piz_ae2b_int(i,k,inu) = piz_aeri_2bands(inu,spinsol)
@@ -676,5 +756,6 @@
 
                  tau_ae(i,k,aerosol_name(m),inu) = &
-                    mass_temp(i,k,7+ spinsol)*1000.*zdp1(i,k)*delt*tau_ae2b_int(i,k,inu)*fac
+                      mass_temp(i,k,naero_soluble+ spinsol)*1000.*zdp1(i,k)* &
+                                               delt*tau_ae2b_int(i,k,inu)*fac
               ENDIF
 
@@ -682,5 +763,4 @@
 
               cg_ae(i,k,aerosol_name(m),inu)= cg_ae2b_int(i,k,inu)
-
 
            ENDDO    ! nbands : boucle sur les bandes spectrale
@@ -726,5 +806,4 @@
                       (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
 
-
               ELSEIF (mrfspecies .EQ. 3) THEN             ! = natural aerosol NAT
 	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)+ &
@@ -767,5 +846,4 @@
                       (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
 
-
               ELSEIF (mrfspecies .EQ. 4) THEN             ! = BC
 	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu)
@@ -778,4 +856,5 @@
                       +tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu))/ &
                       (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
+
               ELSEIF (mrfspecies .EQ. 5) THEN             ! = SO4
 	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu)
@@ -799,4 +878,5 @@
                       +tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu))/ &
                       (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
+
               ELSEIF (mrfspecies .EQ. 7) THEN             ! = DUST
 	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_CIDUSTM,inu)
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90	(revision 1237)
@@ -1,3 +1,5 @@
-
+!
+! $Id$
+!
 
 SUBROUTINE AEROPT_5WV(&
@@ -9,5 +11,5 @@
   USE DIMPHY
   USE aero_mod
-
+  USE mod_phys_lmdz_para, ONLY : mpi_rank
   !
   !    Yves Balkanski le 12 avril 2006
@@ -72,5 +74,5 @@
   LOGICAL :: soluble
   
-  INTEGER :: i, k, m
+  INTEGER :: i, k, ierr, m
   INTEGER :: spsol, spinsol, spss, la
   INTEGER :: RH_num
@@ -81,6 +83,7 @@
   INTEGER, PARAMETER :: la865 = 5
   INTEGER, PARAMETER :: nbre_RH=12
-  INTEGER, PARAMETER :: nbsol_compaer=3   !  1- BC soluble; 2- POM soluble; 3- SO4.
-  INTEGER, PARAMETER :: nbinsol_compaer=3 !  1- Dust; 2- BC insoluble; 3- POM insoluble
+  INTEGER, PARAMETER :: naero_soluble=7   !  1- BC soluble; 2- POM soluble; 3- SO4 acc.
+                                          ! 4- SO4 coarse; 5 seasalt super-C; 6 seasalt coarse; 7 seasalt acc.
+  INTEGER, PARAMETER :: naero_insoluble=3 !  1- Dust; 2- BC insoluble; 3- POM insoluble
   INTEGER, PARAMETER :: nb_level = 19 ! number of vertical levels
   LOGICAL, SAVE :: firstcall=.TRUE.
@@ -98,5 +101,6 @@
 
   ! Coefficient optiques interpole sur le nombre de niveau du modele
-  REAL, DIMENSION(klev) :: A1_ASSSM, A2_ASSSM, A3_ASSSM,&
+  REAL, ALLOCATABLE,  DIMENSION(:), SAVE :: &
+          A1_ASSSM, A2_ASSSM, A3_ASSSM,&
           B1_ASSSM, B2_ASSSM, C1_ASSSM, C2_ASSSM,&
           A1_CSSSM, A2_CSSSM, A3_CSSSM,&
@@ -124,13 +128,13 @@
 
   
-  REAL :: alpha_aers_5wv(nbre_RH,las,nbsol_compaer)   ! ext. coeff. Soluble comp. units *** m2/g 
-                                                      ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4.
-  REAL :: alpha_aeri_5wv(las,nbinsol_compaer)         ! ext. coeff. Insoluble comp. 1- Dust: 2- BC; 3- POM
-  REAL :: cg_aers_5wv(nbre_RH,las,nbsol_compaer)      ! Asym. param. soluble comp. 
-                                                      ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4.
-  REAL :: cg_aeri_5wv(las,nbinsol_compaer)            ! Asym. param. insoluble comp. 1- Dust: 2- BC; 3- POM
-  REAL :: piz_aers_5wv(nbre_RH,las,nbsol_compaer)   
-                                                      ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4.
-  REAL :: piz_aeri_5wv(las,nbinsol_compaer)           ! Insoluble comp. 1- Dust: 2- BC; 3- POM
+  REAL :: alpha_aers_5wv(nbre_RH,las,naero_soluble)   ! ext. coeff. Soluble comp. units *** m2/g 
+   !  1- BC soluble; 2- POM soluble; 3- SO4 acc.; 4- SO4 coarse; 5 seasalt super-C; 6 seasalt coarse; 7 seasalt acc.
+  REAL :: alpha_aeri_5wv(las,naero_insoluble)         ! ext. coeff. Insoluble comp. 1- Dust: 2- BC; 3- POM
+  REAL :: cg_aers_5wv(nbre_RH,las,naero_soluble)      ! Asym. param. soluble comp. 
+   !  1- BC soluble; 2- POM soluble; 3- SO4 acc.; 4- SO4 coarse; 5 seasalt super-C; 6 seasalt coarse; 7 seasalt acc.
+  REAL :: cg_aeri_5wv(las,naero_insoluble)            ! Asym. param. insoluble comp. 1- Dust: 2- BC; 3- POM
+  REAL :: piz_aers_5wv(nbre_RH,las,naero_soluble)   
+   !  1- BC soluble; 2- POM soluble; 3- SO4 acc.; 4- SO4 coarse; 5 seasalt super-C; 6 seasalt coarse; 7 seasalt acc.
+  REAL :: piz_aeri_5wv(las,naero_insoluble)           ! Insoluble comp. 1- Dust: 2- BC; 3- POM
 
   REAL, DIMENSION(klon,klev,naero_spc) :: mass_temp
@@ -284,82 +288,169 @@
   !                                   le 12 AVRIL 2006
   ! 
-  DATA alpha_aers_5wv/ &
-     ! bc soluble
-     7.930,7.930,7.930,7.930,7.930,7.930,     &
-     7.930,7.930,10.893,12.618,14.550,16.613, &
-     7.658,7.658,7.658,7.658,7.658,7.658,     &
-     7.658,7.658,10.351,11.879,13.642,15.510, &
-     7.195,7.195,7.195,7.195,7.195,7.195,     &
-     7.195,7.195,9.551,10.847,12.381,13.994,  &
-     6.736,6.736,6.736,6.736,6.736,6.736,     &
-     6.736,6.736,8.818,9.938,11.283,12.687,   &
-     6.277,6.277,6.277,6.277,6.277,6.277,     &
-     6.277,6.277,8.123,9.094,10.275,11.501,   &
-     ! pom soluble
-     6.676,6.676,6.676,6.676,6.710,6.934,   &
-     7.141,7.569,8.034,8.529,9.456,10.511,  &
-     5.109,5.109,5.109,5.109,5.189,5.535,   &
-     5.960,6.852,8.008,9.712,12.897,19.676, &
-     3.718,3.718,3.718,3.718,3.779,4.042,   &
-     4.364,5.052,5.956,7.314,9.896,15.688,  &
-     2.849,2.849,2.849,2.849,2.897,3.107,   &
-     3.365,3.916,4.649,5.760,7.900,12.863,  &
-     2.229,2.229,2.229,2.229,2.268,2.437,   &
-     2.645,3.095,3.692,4.608,6.391,10.633,  &
-     ! Sulfate
-     5.751,6.215,6.690,7.024,7.599,8.195,      &
-     9.156,10.355,12.660,14.823,18.908,24.508, &
-     4.320,4.675,5.052,5.375,5.787,6.274,      &
-     7.066,8.083,10.088,12.003,15.697,21.133,  &
-     3.079,3.351,3.639,3.886,4.205,4.584,      &
-     5.206,6.019,7.648,9.234,12.391,17.220,    &
-     2.336,2.552,2.781,2.979,3.236,3.540,      &
-     4.046,4.711,6.056,7.388,10.093,14.313,    &
-     1.777,1.949,2.134,2.292,2.503,2.751,      &
-     3.166,3.712,4.828,5.949,8.264,11.922/
+ DATA alpha_aers_5wv/ &
+                                ! bc soluble
+       7.930,7.930,7.930,7.930,7.930,7.930,     &
+       7.930,7.930,10.893,12.618,14.550,16.613, &
+       7.658,7.658,7.658,7.658,7.658,7.658,     &
+       7.658,7.658,10.351,11.879,13.642,15.510, &
+       7.195,7.195,7.195,7.195,7.195,7.195,     &
+       7.195,7.195,9.551,10.847,12.381,13.994,  &
+       6.736,6.736,6.736,6.736,6.736,6.736,     &
+       6.736,6.736,8.818,9.938,11.283,12.687,   &
+       6.277,6.277,6.277,6.277,6.277,6.277,     &
+       6.277,6.277,8.123,9.094,10.275,11.501,   &
+                                ! pom soluble
+       6.676,6.676,6.676,6.676,6.710,6.934,   &
+       7.141,7.569,8.034,8.529,9.456,10.511,  &
+       5.109,5.109,5.109,5.109,5.189,5.535,   &
+       5.960,6.852,8.008,9.712,12.897,19.676, &
+       3.718,3.718,3.718,3.718,3.779,4.042,   &
+       4.364,5.052,5.956,7.314,9.896,15.688,  &
+       2.849,2.849,2.849,2.849,2.897,3.107,   &
+       3.365,3.916,4.649,5.760,7.900,12.863,  &
+       2.229,2.229,2.229,2.229,2.268,2.437,   &
+       2.645,3.095,3.692,4.608,6.391,10.633,  &
+                                ! Sulfate (Accumulation)
+       5.751,6.215,6.690,7.024,7.599,8.195,      &
+       9.156,10.355,12.660,14.823,18.908,24.508, &
+       4.320,4.675,5.052,5.375,5.787,6.274,      &
+       7.066,8.083,10.088,12.003,15.697,21.133,  &
+       3.079,3.351,3.639,3.886,4.205,4.584,      &
+       5.206,6.019,7.648,9.234,12.391,17.220,    &
+       2.336,2.552,2.781,2.979,3.236,3.540,      &
+       4.046,4.711,6.056,7.388,10.093,14.313,    &
+       1.777,1.949,2.134,2.292,2.503,2.751,      &
+       3.166,3.712,4.828,5.949,8.264,11.922,     &
+                                ! Sulfate (Coarse)
+       5.751,6.215,6.690,7.024,7.599,8.195,      &
+       9.156,10.355,12.660,14.823,18.908,24.508, &
+       4.320,4.675,5.052,5.375,5.787,6.274,      &
+       7.066,8.083,10.088,12.003,15.697,21.133,  &
+       3.079,3.351,3.639,3.886,4.205,4.584,      &
+       5.206,6.019,7.648,9.234,12.391,17.220,    &
+       2.336,2.552,2.781,2.979,3.236,3.540,      &
+       4.046,4.711,6.056,7.388,10.093,14.313,    &
+       1.777,1.949,2.134,2.292,2.503,2.751,      &
+       3.166,3.712,4.828,5.949,8.264,11.922,     &
+                                ! Seasalt soluble super_coarse (computed below for 550nm)
+       0.50,0.90,1.05,1.21,1.40,2.41, & 
+       2.66,3.11,3.88,4.52,5.69,8.84, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+     0.52,0.93,1.08,1.24,1.43,2.47, & 
+     2.73,3.20,3.99,4.64,5.84,9.04, & 
+     0.52,0.93,1.09,1.25,1.44,2.50, & 
+     2.76,3.23,4.03,4.68,5.89,9.14, & 
+     0.52,0.94,1.09,1.26,1.45,2.51, & 
+     2.78,3.25,4.06,4.72,5.94,9.22, & 
+                                ! seasalt soluble coarse (computed below for 550nm)
+       0.50,0.90,1.05,1.21,1.40,2.41, & 
+       2.66,3.11,3.88,4.52,5.69,8.84, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+     0.52,0.93,1.08,1.24,1.43,2.47, & 
+     2.73,3.20,3.99,4.64,5.84,9.04, & 
+     0.52,0.93,1.09,1.25,1.44,2.50, & 
+     2.76,3.23,4.03,4.68,5.89,9.14, & 
+     0.52,0.94,1.09,1.26,1.45,2.51, & 
+     2.78,3.25,4.06,4.72,5.94,9.22, & 
+                                ! seasalt soluble accumulation (computed below for 550nm)
+     4.28, 7.17, 8.44, 9.85,11.60,22.44,  & 
+     25.34,30.54,39.38,46.52,59.33,91.77, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+     2.48, 4.22, 5.02, 5.94, 7.11,15.29,  & 
+     17.70,22.31,30.73,38.06,52.15,90.59, & 
+     1.90, 3.29, 3.94, 4.69, 5.65, 12.58, & 
+     14.68,18.77,26.41,33.25,46.77,85.50, & 
+     1.47, 2.59, 3.12, 3.74, 4.54, 10.42, & 
+     12.24,15.82,22.66,28.91,41.54,79.33/
 
   DATA alpha_aeri_5wv/ &
-     ! dust insoluble
-     0.759, 0.770, 0.775, 0.775, 0.772, &
-     !!jb bc insoluble
-     11.536,10.033, 8.422, 7.234, 6.270, &
-     ! pom insoluble
-     5.042, 3.101, 1.890, 1.294, 0.934/
-
-  DATA cg_aers_5wv/ & 
-     ! bc soluble
-     .651, .651, .651, .651, .651, .651, &
-     .651, .651, .738, .764, .785, .800, &
-     .597, .597, .597, .597, .597, .597, &
-     .597, .597, .695, .725, .751, .770, &
-     .543, .543, .543, .543, .543, .543, &
-     .543, .543, .650, .684, .714, .736, & 
-     .504, .504, .504, .504, .504, .504, &
-     .504, .504, .614, .651, .683, .708, & 
-     .469, .469, .469, .469, .469, .469, &
-     .469, .469, .582, .620, .655, .681, &
-     ! pom soluble
-     .679, .679, .679, .679, .683, .691, &
-     .703, .720, .736, .751, .766, .784, &
-     .656, .656, .656, .656, .659, .669, &
-     .681, .699, .717, .735, .750, .779, & 
-     .623, .623, .623, .623, .627, .637, &
-     .649, .668, .688, .709, .734, .762, &
-     .592, .592, .592, .592, .595, .605, &
-     .618, .639, .660, .682, .711, .743, &
-     .561, .561, .561, .561, .565, .575, &
-     .588, .609, .632, .656, .688, .724, &
-     ! sulfate
-     .671, .684, .697, .704, .714, .723, &
-     .734, .746, .762, .771, .781, .789, &
-     .653, .666, .678, .687, .697, .707, &
-     .719, .732, .751, .762, .775, .789, &
-     .622, .635, .648, .657, .667, .678, &
-     .691, .705, .728, .741, .758, .777, &
-     .591, .604, .617, .627, .638, .650, &
-     .664, .679, .704, .719, .739, .761, &
-     .560, .574, .587, .597, .609, .621, & 
-     .637, .653, .680, .697, .719, .745/
-  !
+                                ! dust insoluble
+       0.759, 0.770, 0.775, 0.775, 0.772, &
+                                !!jb bc insoluble
+       11.536,10.033, 8.422, 7.234, 6.270, &
+                                ! pom insoluble
+       5.042, 3.101, 1.890, 1.294, 0.934/
+  !
+ DATA cg_aers_5wv/ & 
+                                ! bc soluble
+       .651, .651, .651, .651, .651, .651, &
+       .651, .651, .738, .764, .785, .800, &
+       .597, .597, .597, .597, .597, .597, &
+       .597, .597, .695, .725, .751, .770, &
+       .543, .543, .543, .543, .543, .543, &
+       .543, .543, .650, .684, .714, .736, & 
+       .504, .504, .504, .504, .504, .504, &
+       .504, .504, .614, .651, .683, .708, & 
+       .469, .469, .469, .469, .469, .469, &
+       .469, .469, .582, .620, .655, .681, &
+                                ! pom soluble
+       .679, .679, .679, .679, .683, .691, &
+       .703, .720, .736, .751, .766, .784, &
+       .656, .656, .656, .656, .659, .669, &
+       .681, .699, .717, .735, .750, .779, & 
+       .623, .623, .623, .623, .627, .637, &
+       .649, .668, .688, .709, .734, .762, &
+       .592, .592, .592, .592, .595, .605, &
+       .618, .639, .660, .682, .711, .743, &
+       .561, .561, .561, .561, .565, .575, &
+       .588, .609, .632, .656, .688, .724, &
+                                ! Accumulation sulfate
+       .671, .684, .697, .704, .714, .723, &
+       .734, .746, .762, .771, .781, .789, &
+       .653, .666, .678, .687, .697, .707, &
+       .719, .732, .751, .762, .775, .789, &
+       .622, .635, .648, .657, .667, .678, &
+       .691, .705, .728, .741, .758, .777, &
+       .591, .604, .617, .627, .638, .650, &
+       .664, .679, .704, .719, .739, .761, &
+       .560, .574, .587, .597, .609, .621, & 
+       .637, .653, .680, .697, .719, .745, &
+                                ! Coarse sulfate
+       .671, .684, .697, .704, .714, .723, &
+       .734, .746, .762, .771, .781, .789, &
+       .653, .666, .678, .687, .697, .707, &
+       .719, .732, .751, .762, .775, .789, &
+       .622, .635, .648, .657, .667, .678, &
+       .691, .705, .728, .741, .758, .777, &
+       .591, .604, .617, .627, .638, .650, &
+       .664, .679, .704, .719, .739, .761, &
+       .560, .574, .587, .597, .609, .621, & 
+       .637, .653, .680, .697, .719, .745, &
+                                ! For super coarse seasalt (computed below for 550nm!)
+      0.730,0.753,0.760,0.766,0.772,0.793, & 
+      0.797,0.802,0.809,0.813,0.820,0.830, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+     0.721,0.744,0.750,0.756,0.762,0.784, & 
+     0.787,0.793,0.800,0.804,0.811,0.822, & 
+     0.717,0.741,0.747,0.753,0.759,0.780, & 
+     0.784,0.789,0.795,0.800,0.806,0.817, & 
+     0.715,0.739,0.745,0.751,0.757,0.777, &  
+     0.781,0.786,0.793,0.797,0.803,0.814, & 
+                                ! For coarse-soluble seasalt (computed below for 550nm!)
+     0.730,0.753,0.760,0.766,0.772,0.793, & 
+     0.797,0.802,0.809,0.813,0.820,0.830, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+     0.721,0.744,0.750,0.756,0.762,0.784, & 
+     0.787,0.793,0.800,0.804,0.811,0.822, & 
+     0.717,0.741,0.747,0.753,0.759,0.780, & 
+     0.784,0.789,0.795,0.800,0.806,0.817, & 
+     0.715,0.739,0.745,0.751,0.757,0.777, &  
+     0.781,0.786,0.793,0.797,0.803,0.814, & 
+                                ! accumulation-seasalt soluble (computed below for 550nm!) 
+     0.698,0.722,0.729,0.736,0.743,0.765, & 
+     0.768,0.773,0.777,0.779,0.781,0.779, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+       0.000,0.000,0.000,0.000,0.000,0.000, & 
+     0.658,0.691,0.701,0.710,0.720,0.756, & 
+     0.763,0.771,0.782,0.788,0.795,0.801, & 
+     0.632,0.668,0.679,0.690,0.701,0.743, & 
+     0.750,0.762,0.775,0.783,0.792,0.804, & 
+     0.605,0.644,0.656,0.669,0.681,0.729, & 
+     0.737,0.750,0.765,0.775,0.787,0.803/
 
   DATA cg_aeri_5wv/&
@@ -372,37 +463,81 @@
   !
   DATA piz_aers_5wv/&
-     ! bc soluble
-     .445, .445, .445, .445, .445, .445, &
-     .445, .445, .470, .487, .508, .531, &
-     .442, .442, .442, .442, .442, .442, &
-     .442, .442, .462, .481, .506, .533, &
-     .427, .427, .427, .427, .427, .427, &
-     .427, .427, .449, .470, .497, .526, &
-     .413, .413, .413, .413, .413, .413, &
-     .413, .413, .437, .458, .486, .516, &
-     .399, .399, .399, .399, .399, .399, &
-     .399, .399, .423, .445, .473, .506, &
-     ! pom soluble
-     .975, .975, .975, .975, .975, .977, &
-     .979, .982, .984, .987, .990, .994, &
-     .972, .972, .972, .972, .973, .974, &
-     .977, .980, .983, .986, .989, .993, &
-     .963, .963, .963, .963, .964, .966, &
-     .969, .974, .977, .982, .986, .991, &
-     .955, .955, .955, .955, .955, .958, &
-     .962, .967, .972, .977, .983, .989, &
-     .944, .944, .944, .944, .944, .948, &
-     .952, .959, .962, .972, .979, .987, &
-     ! sulfate
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000, &
-     1.000,1.000,1.000,1.000,1.000,1.000/
+                                ! bc soluble
+       .445, .445, .445, .445, .445, .445, &
+       .445, .445, .470, .487, .508, .531, &
+       .442, .442, .442, .442, .442, .442, &
+       .442, .442, .462, .481, .506, .533, &
+       .427, .427, .427, .427, .427, .427, &
+       .427, .427, .449, .470, .497, .526, &
+       .413, .413, .413, .413, .413, .413, &
+       .413, .413, .437, .458, .486, .516, &
+       .399, .399, .399, .399, .399, .399, &
+       .399, .399, .423, .445, .473, .506, &
+                                ! pom soluble
+       .975, .975, .975, .975, .975, .977, &
+       .979, .982, .984, .987, .990, .994, &
+       .972, .972, .972, .972, .973, .974, &
+       .977, .980, .983, .986, .989, .993, &
+       .963, .963, .963, .963, .964, .966, &
+       .969, .974, .977, .982, .986, .991, &
+       .955, .955, .955, .955, .955, .958, &
+       .962, .967, .972, .977, .983, .989, &
+       .944, .944, .944, .944, .944, .948, &
+       .952, .959, .962, .972, .979, .987, &
+                                ! sulfate soluble accumulation
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+                                ! sulfate soluble coarse
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+                                ! seasalt super coarse (computed below for 550nm)
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+                                ! seasalt coarse (computed below for 550nm)
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+                                ! seasalt soluble accumulation (computed below for 550nm)
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000, & 
+       1.000,1.000,1.000,1.000,1.000,1.000/
   !
   DATA piz_aeri_5wv/&
@@ -417,4 +552,15 @@
   IF (firstcall) THEN
      firstcall=.FALSE.
+! Allocation
+    IF (.NOT. ALLOCATED(A1_ASSSM)) THEN
+        ALLOCATE(A1_ASSSM(klev),A2_ASSSM(klev), A3_ASSSM(klev),&
+          B1_ASSSM(klev), B2_ASSSM(klev), C1_ASSSM(klev), C2_ASSSM(klev),&
+          A1_CSSSM(klev), A2_CSSSM(klev), A3_CSSSM(klev),&
+          B1_CSSSM(klev), B2_CSSSM(klev), C1_CSSSM(klev), C2_CSSSM(klev),&
+          A1_SSSSM(klev), A2_SSSSM(klev), A3_SSSSM(klev),&
+          B1_SSSSM(klev), B2_SSSSM(klev), C1_SSSSM(klev), C2_SSSSM(klev), stat=ierr)
+        IF (ierr /= 0) CALL abort_gcm('aeropt_5mw', 'pb in allocation 1',1)
+     END IF
+
 !Accumulation mode
      CALL pres2lev(A1_ASSSM_19, A1_ASSSM, nb_level, klev, presnivs_19, presnivs, 1, 1, .FALSE.)
@@ -456,10 +602,8 @@
   DO k=1, klev
     DO i=1, klon
-      IF (t_seri(i,k).EQ.0) stop 'stop aeropt_5wv T '
-      IF (pplay(i,k).EQ.0) stop  'stop aeropt_5wv p '
-
+!      IF (t_seri(i,k).EQ.0) stop 'stop aeropt_5wv T '
+!      IF (pplay(i,k).EQ.0) stop  'stop aeropt_5wv p '
       zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
       mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9
-
     ENDDO
   ENDDO
@@ -519,6 +663,5 @@
   zdp1=pdel(:,:)/(gravit*delt) 
   
-  DO m=1,nb_aer   ! tau is only computed for each mass
-    
+  DO m=1,nb_aer   ! tau is only computed for each mass    
     fac=1.0
     IF (aerosol_name(m).EQ.id_ASBCM) THEN
@@ -530,7 +673,12 @@
         spsol=2 
         spss=0
-    ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR.  (aerosol_name(m).EQ.id_CSSO4M)) THEN
+    ELSEIF (aerosol_name(m).EQ.id_ASSO4M) THEN
         soluble=.TRUE.
         spsol=3
+        spss=0
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+    ELSEIF (aerosol_name(m).EQ.id_CSSO4M) THEN
+        soluble=.TRUE.
+        spsol=4
         spss=0
         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
@@ -581,5 +729,4 @@
           IF (soluble) THEN
 
-
               IF((la.EQ.2).AND.(spss.NE.0)) THEN !la=2 corresponds to 550 nm
                   H=rh/100
@@ -620,5 +767,8 @@
                  mass_temp(i,k,spsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac
 
-          ELSE
+               tausum(i,la,spsol)=tausum(i,la,spsol)+tau3d(i,k)
+           tau(i,k,la,spsol)=tau3d(i,k)
+
+          ELSE                                           ! Case insoluble aerosol
               tau_ae5wv_int(i,k,la) = alpha_aeri_5wv(la,spinsol)
               piz_ae5wv_int(i,k,la) = piz_aeri_5wv(la,spinsol)
@@ -626,32 +776,14 @@
 
               tau3d(i,k) = &
-                 mass_temp(i,k,7+spinsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac
+                 mass_temp(i,k,naero_soluble+spinsol)*1000.*zdp1(i,k)* &
+                      tau_ae5wv_int(i,k,la)*delt*fac
+
+               tausum(i,la,naero_soluble+spinsol)=tausum(i,la,spsol)+tau3d(i,k)
+           tau(i,k,la,naero_soluble+spinsol)=tau3d(i,k)
+
           ENDIF
-          
-          
+
         ENDDO     ! Boucle sur les points géographiques (grille horizontale)
       ENDDO     ! Boucle sur les niveaux verticaux
-
-      IF (soluble) THEN
-
-          tau(:,:,la,spsol)=tau3d(:,:)
-      
-          DO k=1, KLEV
-            DO i=1,KLON
-              tausum(i,la,spsol)=tausum(i,la,spsol)+tau3d(i,k)
-            ENDDO
-          ENDDO
-      ELSE
-          tau(:,:,la,spsol)=tau3d(:,:)
-      
-          DO k=1, KLEV
-            DO i=1,KLON
-              tausum(i,la,5+spinsol)=tausum(i,la,5+spinsol)+tau3d(i,k)
-            ENDDO
-          ENDDO
-      ENDIF
-
-
-
     ENDDO   ! boucle sur les longueurs d'onde
   ENDDO     ! Boucle  sur les masses de traceurs
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/hines_gwd.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/hines_gwd.F	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/hines_gwd.F	(revision 1237)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
       SUBROUTINE HINES_GWD(NLON,NLEV,DTIME,paphm1x, papm1x,
      I      rlat,tx,ux,vx,
@@ -1666,23 +1669,23 @@
 C  the variances.
 C
-      DO 80 N = 1,NAZ
-        DO 70 I = IL1,IL2
-          IF (I_ALPHA(I,N).LT.0.)  THEN
-            WRITE (6,*) 
-            WRITE (6,*) '******************************'
-            WRITE (6,*) 'Hines integral I_ALPHA < 0 '
-            WRITE (6,*) '  longitude I=',I
-            WRITE (6,*) '  azimuth   N=',N
-            WRITE (6,*) '  level   LEV=',LEV
-            WRITE (6,*) '  I_ALPHA =',I_ALPHA(I,N)
-            WRITE (6,*) '  V_ALPHA =',V_ALPHA(I,LEV,N)
-            WRITE (6,*) '  M_ALPHA =',M_ALPHA(I,LEV,N)
-            WRITE (6,*) '  Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I)
-            WRITE (6,*) '  QM      =',V_ALPHA(I,LEV,N) / BVFB(I) 
-     ^                                * M_ALPHA(I,LEV,N)
-            WRITE (6,*) '******************************'
-          END IF
- 70     CONTINUE
- 80   CONTINUE
+c      DO 80 N = 1,NAZ
+c        DO 70 I = IL1,IL2
+c          IF (I_ALPHA(I,N).LT.0.)  THEN
+c            WRITE (6,*) 
+c            WRITE (6,*) '******************************'
+c            WRITE (6,*) 'Hines integral I_ALPHA < 0 '
+c            WRITE (6,*) '  longitude I=',I
+c            WRITE (6,*) '  azimuth   N=',N
+c            WRITE (6,*) '  level   LEV=',LEV
+c            WRITE (6,*) '  I_ALPHA =',I_ALPHA(I,N)
+c            WRITE (6,*) '  V_ALPHA =',V_ALPHA(I,LEV,N)
+c            WRITE (6,*) '  M_ALPHA =',M_ALPHA(I,LEV,N)
+c            WRITE (6,*) '  Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I)
+c            WRITE (6,*) '  QM      =',V_ALPHA(I,LEV,N) / BVFB(I) 
+c     ^                                * M_ALPHA(I,LEV,N)
+c            WRITE (6,*) '******************************'
+c          END IF
+c 70     CONTINUE
+c 80   CONTINUE
 C
       RETURN
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F	(revision 1237)
@@ -164,5 +164,5 @@
                                 !
                   cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
-     &               log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
+     &             log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
                   cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
                ENDDO
@@ -178,5 +178,5 @@
      &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k)))  
      &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
-                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.) 
+                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) *1.e6, 5.) 
                ENDDO            
             ENDDO
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1237)
@@ -2671,6 +2671,6 @@
          IF (.NOT. aerosol_couple)
      &        CALL readaerosol_optic(
-     &        debut, new_aod, flag_aerosol, jD_cur-jD_ref, pdtphys,
-     &        pplay, paprs, t_seri, rhcl, presnivs, 
+     &        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
+     &        pdtphys, pplay, paprs, t_seri, rhcl, presnivs, 
      &        mass_solu_aero, mass_solu_aero_pi, 
      &        tau_aero, piz_aero, cg_aero, 
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90	(revision 1237)
@@ -1,5 +1,5 @@
 ! $Id$
 !
-SUBROUTINE readaerosol_interp(id_aero, r_day, first, pplay, paprs, t_seri, mass_out, pi_mass_out)
+SUBROUTINE readaerosol_interp(id_aero, itap, pdtphys, r_day, first, pplay, paprs, t_seri, mass_out, pi_mass_out)
 !
 ! This routine will return the mass concentration at actual day(mass_out) and 
@@ -18,4 +18,5 @@
   USE aero_mod, ONLY : naero_spc, name_aero
   USE write_field_phy
+  USE phys_cal_mod
 
   IMPLICIT NONE
@@ -32,4 +33,6 @@
 !****************************************************************************************
   INTEGER, INTENT(IN)                    :: id_aero! Identity number for the aerosol to treat
+  INTEGER, INTENT(IN)                    :: itap   ! Physic step count
+  REAL, INTENT(IN)                       :: pdtphys! Physic day step
   REAL, INTENT(IN)                       :: r_day  ! Day of integration
   LOGICAL, INTENT(IN)                    :: first  ! First model timestep 
@@ -46,5 +49,5 @@
 !****************************************************************************************
   INTEGER                         :: i, k, ierr
-  INTEGER                         :: iday, iyr
+  INTEGER                         :: iday, iyr, lmt_pas
   INTEGER                         :: im, day1, day2, im2
   INTEGER                         :: pi_klev_src ! Only for testing purpose
@@ -77,4 +80,5 @@
 
   LOGICAL            :: lnewday      ! Indicates if first time step at a new day
+  LOGICAL            :: OLDNEWDAY
   LOGICAL,SAVE       :: vert_interp  ! Indicates if vertical interpolation will be done
   LOGICAL,SAVE       :: debug=.FALSE.! Debugging in this subroutine
@@ -88,5 +92,15 @@
 
 ! Calculation to find if it is a new day
-  iday = INT(r_day) 
+
+  IF(mpi_rank == 0)then
+     PRINT*,'CONTROL PANEL REGARDING TIME STEPING'
+  ENDIF
+
+  ! Use phys_cal_mod
+  !iday= day_cur
+  !iyr = year_cur
+  !im  = mth_cur
+
+  iday = INT(r_day)
   iyr  = iday/360
   iday = iday-iyr*360         ! day of the actual year
@@ -94,6 +108,24 @@
   im   = iday/30 +1           ! the actual month
 
-  ! 0.02 is about 0.5/24, namly less than half an hour
-  lnewday = (r_day-FLOAT(iday) < 0.02)
+  IF(MOD(itap-1,NINT(86400./pdtphys)) == 0)THEN
+     lnewday=.TRUE.
+  ENDIF
+
+  IF(mpi_rank == 0)then
+     ! 0.02 is about 0.5/24, namly less than half an hour
+     OLDNEWDAY = (r_day-FLOAT(iday) < 0.02)
+     ! Once per day, update aerosol fields
+     lmt_pas = NINT(86400./pdtphys)
+     PRINT*,'r_day-FLOAT(iday) =',r_day-FLOAT(iday) 
+     PRINT*,'itap =',itap
+     PRINT*,'pdtphys =',pdtphys
+     PRINT*,'lmt_pas =',lmt_pas
+     PRINT*,'iday =',iday
+     PRINT*,'r_day =',r_day
+     PRINT*,'NINT(86400./pdtphys) =',NINT(86400./pdtphys)
+     PRINT*,'MOD(0,1) =',MOD(0,1)
+     PRINT*,'lnewday =',lnewday
+     PRINT*,'OLDNEWDAY =',OLDNEWDAY
+  ENDIF
 
   IF (.NOT. ALLOCATED(var_day)) THEN
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_optic.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_optic.F90	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_optic.F90	(revision 1237)
@@ -1,6 +1,6 @@
 ! $Id$
 !
-SUBROUTINE readaerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, &
-     pplay, paprs, t_seri, rhcl, presnivs, &
+SUBROUTINE readaerosol_optic(debut, new_aod, flag_aerosol, itap, rjourvrai, &
+     pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
      mass_solu_aero, mass_solu_aero_pi, &
      tau_aero, piz_aero, cg_aero, &
@@ -21,4 +21,5 @@
   LOGICAL, INTENT(IN)                      :: new_aod
   INTEGER, INTENT(IN)                      :: flag_aerosol
+  INTEGER, INTENT(IN)                      :: itap
   REAL, INTENT(IN)                         :: rjourvrai
   REAL, INTENT(IN)                         :: pdtphys
@@ -74,5 +75,5 @@
        flag_aerosol .EQ. 6 ) THEN 
 
-     CALL readaerosol_interp(id_ASSO4M, rjourvrai, debut, pplay, paprs, t_seri, sulfate, sulfate_pi)
+     CALL readaerosol_interp(id_ASSO4M, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sulfate, sulfate_pi)
   ELSE
      sulfate(:,:) = 0. ; sulfate_pi(:,:) = 0.
@@ -84,6 +85,6 @@
 
      ! Get bc aerosol distribution 
-     CALL readaerosol_interp(id_ASBCM, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi )
-     CALL readaerosol_interp(id_AIBCM, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi )
+     CALL readaerosol_interp(id_ASBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi )
+     CALL readaerosol_interp(id_AIBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi )
   ELSE
      bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0.
@@ -96,6 +97,6 @@
        flag_aerosol .EQ. 6 ) THEN
 
-     CALL readaerosol_interp(id_ASPOMM, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi)
-     CALL readaerosol_interp(id_AIPOMM, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi)
+     CALL readaerosol_interp(id_ASPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi)
+     CALL readaerosol_interp(id_AIPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi)
   ELSE
      pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
@@ -108,7 +109,7 @@
       flag_aerosol .EQ. 6 ) THEN 
 
-      CALL readaerosol_interp(id_SSSSM ,rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi) 
-      CALL readaerosol_interp(id_CSSSM ,rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi) 
-      CALL readaerosol_interp(id_ASSSM ,rjourvrai, debut, pplay, paprs, t_seri, ssacu, ssacu_pi) 
+      CALL readaerosol_interp(id_SSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi) 
+      CALL readaerosol_interp(id_CSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi) 
+      CALL readaerosol_interp(id_ASSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, ssacu, ssacu_pi) 
 
   ELSE
@@ -122,5 +123,5 @@
       flag_aerosol .EQ. 6 ) THEN 
 
-      CALL readaerosol_interp(id_CIDUSTM, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi) 
+      CALL readaerosol_interp(id_CIDUSTM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi) 
 
   ELSE
@@ -163,5 +164,5 @@
 
      fractnat_allaer(:,:) = 0.
-     CALL aeropt_2bands( &
+     CALL aeropt_2bands(                 &
           pdel, m_allaer, pdtphys, rhcl, & 
           tau_aero, piz_aero, cg_aero,   &
@@ -170,8 +171,8 @@
      
      ! aeropt_5wv only for validation and diagnostics.
-     CALL aeropt_5wv( &
-          pdel, m_allaer, &
-          pdtphys, rhcl, aerindex, & 
-          flag_aerosol, pplay, t_seri, &
+     CALL aeropt_5wv(                    &
+          pdel, m_allaer,                &
+          pdtphys, rhcl, aerindex,       & 
+          flag_aerosol, pplay, t_seri,   &
           tausum_aero, tau3d_aero, presnivs)
   ELSE
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aeroAR4.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aeroAR4.F90	(revision 1235)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aeroAR4.F90	(revision 1237)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 SUBROUTINE SW_AEROAR4(PSCT, PRMU0, PFRAC, &
      PPMB, PDP, &
@@ -169,4 +172,6 @@
   !$OMP THREADPRIVATE(ZFSDN0_AERO)
 
+!
+  LOGICAL :: AEROSOLFEEDBACK_ACTIVE=.true.
 
   IF(.NOT.initialized) THEN
@@ -207,5 +212,6 @@
         DO JL = 1, KDLON
            ZCLDSW0(JL,JK) = 0.0
-           ZOZ(JL,JK) = POZON(JL,JK) / dobson_u / 1e3 / RG * PDP(JL,JK)
+           ZOZ(JL,JK) = POZON(JL,JK)*46.6968/RG &
+                *PDP(JL,JK)*(101325.0/PPSOL(JL))
         ENDDO
      ENDDO
@@ -232,8 +238,6 @@
      DO JK = 1 , KFLEV+1
         DO JL = 1, KDLON
-           ZFSUP0(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-           ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ZFSUP0_AERO(JL,JK,1) = ZFSUP0(JL,JK) 
-           ZFSDN0_AERO(JL,JK,1) = ZFSDN0(JL,JK)
+           ZFSUP0_AERO(JL,JK,1) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+           ZFSDN0_AERO(JL,JK,1) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
         ENDDO
      ENDDO
@@ -261,15 +265,9 @@
      DO JK = 1 , KFLEV+1
         DO JL = 1, KDLON
-           ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-           ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ZFSUP_AERO(JL,JK,1) = ZFSUP(JL,JK) 
-           ZFSDN_AERO(JL,JK,1) = ZFSDN(JL,JK)
+           ZFSUP_AERO(JL,JK,1) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+           ZFSDN_AERO(JL,JK,1) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
         ENDDO
      ENDDO
 
-     ZFSUP0_AERO(:,:,2) = ZFSUP0_AERO(:,:,1)
-     ZFSDN0_AERO(:,:,2) = ZFSDN0_AERO(:,:,1)
-     ZFSUP_AERO(:,:,2) = ZFSUP_AERO(:,:,1)
-     ZFSDN_AERO(:,:,2) = ZFSDN_AERO(:,:,1)
 
 
@@ -298,10 +296,6 @@
         DO JK = 1 , KFLEV+1
            DO JL = 1, KDLON
-              ZFSUPAD0_AERO(JL,JK) = ZFSUP0(JL,JK) 
-              ZFSDNAD0_AERO(JL,JK) = ZFSDN0(JL,JK) 
-              ZFSUP0(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-              ZFSUP0_AERO(JL,JK,2) = ZFSUP0(JL,JK) 
-              ZFSDN0_AERO(JL,JK,2) = ZFSDN0(JL,JK)
+              ZFSUP0_AERO(JL,JK,2) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL) 
+              ZFSDN0_AERO(JL,JK,2) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) 
            ENDDO
         ENDDO
@@ -329,12 +323,10 @@
         DO JK = 1 , KFLEV+1
            DO JL = 1, KDLON
-              ZFSUPAD_AERO(JL,JK) = ZFSUP(JL,JK) 
-              ZFSDNAD_AERO(JL,JK) = ZFSDN(JL,JK) 
-              ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-              ZFSUP_AERO(JL,JK,2) = ZFSUP(JL,JK) 
-              ZFSDN_AERO(JL,JK,2) = ZFSDN(JL,JK)
+              ZFSUP_AERO(JL,JK,2) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL) 
+              ZFSDN_AERO(JL,JK,2) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) 
            ENDDO
         ENDDO
+
+
 
         !CAS NAT
@@ -358,4 +350,5 @@
              ZFDOWN, ZFUP)
 
+! Natural aerosol fluxes
         DO JK = 1 , KFLEV+1
            DO JL = 1, KDLON
@@ -392,273 +385,4 @@
         ENDDO
 
-        ! clear sky (Yves 01/12/2007)
-        ! CAS BC (4)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP0_AERO(JL,JK,4) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN0_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! cloudy-sky + aerosol dir OB
-        ! CAS BC (4)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP_AERO(JL,JK,4) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! clear sky (Yves 13/12/2007)
-        ! CAS SO4 (5)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP0_AERO(JL,JK,5) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN0_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! cloudy-sky + aerosol dir OB
-        ! CAS SO4 (5)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP_AERO(JL,JK,5) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! clear sky (Yves 13/12/2007)
-        ! CAS POM (6)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP0_AERO(JL,JK,6) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN0_AERO(JL,JK,6) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! cloudy-sky + aerosol dir OB
-        ! CAS POM (6)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP_AERO(JL,JK,6) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN_AERO(JL,JK,6) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! clear sky (Yves 13/12/2007)
-        ! CAS DUST (7)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP0_AERO(JL,JK,7) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN0_AERO(JL,JK,7) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! cloudy-sky + aerosol dir OB
-        ! CAS DUST (7)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP_AERO(JL,JK,7) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN_AERO(JL,JK,7) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! clear sky (Yves 13/12/2007)
-        ! CAS Seasalt SS (8)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP0_AERO(JL,JK,8) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN0_AERO(JL,JK,8) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
-
-        ! cloudy-sky + aerosol dir OB
-        ! CAS Seasalt SS (8)
-        flag_aer=1.0
-        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
-             PRMU0,PFRAC,PTAVE,PWV,&
-             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
-        INU = 1
-        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
-             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             ZFD, ZFU)
-        INU = 2
-        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
-             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
-             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
-             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
-             PWV, PQS,&
-             ZFDOWN, ZFUP)
-
-        DO JK = 1 , KFLEV+1
-           DO JL = 1, KDLON
-              ZFSUP_AERO(JL,JK,8) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN_AERO(JL,JK,8) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-           ENDDO
-        ENDDO
 
      ENDIF ! ok_ade
@@ -686,10 +410,6 @@
         DO JK = 1 , KFLEV+1
            DO JL = 1, KDLON
-              ZFSUPAI_AERO(JL,JK) = ZFSUP(JL,JK) 
-              ZFSDNAI_AERO(JL,JK) = ZFSDN(JL,JK)          
-              ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
-              ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
-              ZFSUP_AERO(JL,JK,2) = ZFSUP(JL,JK) 
-              ZFSDN_AERO(JL,JK,2) = ZFSDN(JL,JK)
+              ZFSUP_AERO(JL,JK,4) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) 
            ENDDO
         ENDDO
@@ -699,4 +419,44 @@
   ENDIF
   itapsw = itapsw + 1
+
+
+  IF ( ok_ade .and. ok_aie  ) THEN
+    ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
+    ZFSDN(:,:) =    ZFSDN_AERO(:,:,4)
+    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,2)
+    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,2)
+  ENDIF
+  IF ( ok_ade .and. (.not. ok_aie) )  THEN
+    ZFSUP(:,:) =    ZFSUP_AERO(:,:,2)
+    ZFSDN(:,:) =    ZFSDN_AERO(:,:,2)
+    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,2)
+    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,2)
+  ENDIF
+!MS the following combination would include the direct aerosol effect in cloud regions
+!   because it takes the total aerosol effect 
+  IF ( (.not. ok_ade) .and. ok_aie  )  THEN
+    ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
+    ZFSDN(:,:) =    ZFSDN_AERO(:,:,4)
+    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,1)
+    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,1)
+  ENDIF
+  IF ((.not. ok_ade) .and. (.not. ok_aie)) THEN
+    ZFSUP(:,:) =    ZFSUP_AERO(:,:,1)
+    ZFSDN(:,:) =    ZFSDN_AERO(:,:,1)
+    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,1)
+    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,1)
+  ENDIF
+
+! MS the following allows to compute the forcing diagostics without
+! letting the aerosol forcing act on the meteorology
+! assuming that the no-aerosol case creates the reference meteorological state
+! for the natural aerosol state use: *_AERO(:,:3)
+  IF  (.not. AEROSOLFEEDBACK_ACTIVE) THEN
+    ZFSUP(:,:) =    ZFSUP_AERO(:,:,1)
+    ZFSDN(:,:) =    ZFSDN_AERO(:,:,1)
+    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,1)
+    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,1)
+  ENDIF
+  
 
   DO k = 1, KFLEV
@@ -704,13 +464,12 @@
      DO i = 1, KDLON
 
-        PHEAT(i,k) = -(ZFSUP_AERO(i,kpl1,2)-ZFSUP_AERO(i,k,2)) &
-             -(ZFSDN_AERO(i,k,2)-ZFSDN_AERO(i,kpl1,2))
+        PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k))-(ZFSDN(i,k)-ZFSDN(i,kpl1))
         PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
-        PHEAT0(i,k) = -(ZFSUP0_AERO(i,kpl1,2)-ZFSUP0_AERO(i,k,2)) &
-             -(ZFSDN0_AERO(i,k,2)-ZFSDN0_AERO(i,kpl1,2))
+        PHEAT0(i,k) = -(ZFSUP0(i,kpl1)-ZFSUP0(i,k))-(ZFSDN0(i,k)-ZFSDN0(i,kpl1))
         PHEAT0(i,k) = PHEAT0(i,k) * RDAY*RG/RCPD / PDP(i,k)
 
      ENDDO
   ENDDO
+
   DO i = 1, KDLON
      PALBPLA(i) = ZFSUP(i,KFLEV+1)/(ZFSDN(i,KFLEV+1)+1.0e-20)
@@ -722,23 +481,26 @@
      PTOPSW(i) = ZFSDN(i,KFLEV+1) - ZFSUP(i,KFLEV+1)
 
-     PSOLSW0AERO(i,:) = ZFSDN0_AERO(i,1,:) - ZFSUP0_AERO(i,1,:)
-     PTOPSW0AERO(i,:) = &
-          ZFSDN0_AERO(i,KFLEV+1,:) - ZFSUP0_AERO(i,KFLEV+1,:)
-
-     PSOLSWAERO(i,:) = ZFSDN_AERO(i,1,:) - ZFSUP_AERO(i,1,:)
-     PTOPSWAERO(i,:) = &
-          ZFSDN_AERO(i,KFLEV+1,:) - ZFSUP_AERO(i,KFLEV+1,:)
-
-     PSOLSWADAERO(i) = ZFSDNAD_AERO(i,1) - ZFSUPAD_AERO(i,1)
-     PTOPSWADAERO(i) = &
-          ZFSDNAD_AERO(i,KFLEV+1) - ZFSUPAD_AERO(i,KFLEV+1)
-
-     PSOLSWAD0AERO(i) = ZFSDNAD0_AERO(i,1) - ZFSUPAD0_AERO(i,1)
-     PTOPSWAD0AERO(i) = &
-          ZFSDNAD0_AERO(i,KFLEV+1) - ZFSUPAD0_AERO(i,KFLEV+1)
-
-     PSOLSWAIAERO(i) = ZFSDNAI_AERO(i,1) - ZFSUPAI_AERO(i,1)
-     PTOPSWAIAERO(i) = &
-          ZFSDNAI_AERO(i,KFLEV+1) - ZFSUPAI_AERO(i,KFLEV+1)
+! MS the following is not output, so commented
+!     PSOLSW0AERO(i,:) = ZFSDN0_AERO(i,1,:) - ZFSUP0_AERO(i,1,:)
+!     PTOPSW0AERO(i,:) = &
+!          ZFSDN0_AERO(i,KFLEV+1,:) - ZFSUP0_AERO(i,KFLEV+1,:)
+
+!     PSOLSWAERO(i,:) = ZFSDN_AERO(i,1,:) - ZFSUP_AERO(i,1,:)
+!     PTOPSWAERO(i,:) = &
+!          ZFSDN_AERO(i,KFLEV+1,:) - ZFSUP_AERO(i,KFLEV+1,:)
+
+
+if (ok_ade) then
+     PSOLSWADAERO(i) = (ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))-(ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))
+     PTOPSWADAERO(i) = (ZFSDN_AERO(i,KFLEV+1,2) - ZFSUP_AERO(i,KFLEV+1,2))- (ZFSDN_AERO(i,KFLEV+1,3) - ZFSUP_AERO(i,KFLEV+1,3))
+
+     PSOLSWAD0AERO(i) = (ZFSDN0_AERO(i,1,2) - ZFSUP0_AERO(i,1,2))-(ZFSDN0_AERO(i,1,3) - ZFSUP0_AERO(i,1,3))
+     PTOPSWAD0AERO(i) = (ZFSDN0_AERO(i,KFLEV+1,2) - ZFSUP0_AERO(i,KFLEV+1,2))-(ZFSDN0_AERO(i,KFLEV+1,3) - ZFSUP0_AERO(i,KFLEV+1,3))
+endif
+
+if (ok_aie) then
+     PSOLSWAIAERO(i) = (ZFSDN_AERO(i,1,4) - ZFSUP_AERO(i,1,4))-(ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))
+     PTOPSWAIAERO(i) = (ZFSDN_AERO(i,KFLEV+1,4) - ZFSUP_AERO(i,KFLEV+1,4))-(ZFSDN_AERO(i,KFLEV+1,2) - ZFSUP_AERO(i,KFLEV+1,2))
+endif
 
   ENDDO
