Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90	(revision 1249)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90	(revision 1250)
@@ -51,7 +51,7 @@
   REAL, DIMENSION(klon,klev,naero_tot,nbands) ::  cg_ae
   LOGICAL ::  soluble
-  INTEGER :: i, k, ierr, inu, m, mrfspecies
+  INTEGER :: i, k,n, ierr, inu, m, mrfspecies
   INTEGER :: spsol, spinsol, spss
-  INTEGER :: RH_num
+  INTEGER :: RH_num(klon,klev)
   INTEGER, PARAMETER :: nb_level=19 ! number of vertical levels in DATA
 
@@ -122,8 +122,9 @@
   REAL,PARAMETER :: RH_tab(nbre_RH)=(/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./)
   REAL, PARAMETER :: RH_MAX=95.
-  REAL:: DELTA, rh, H
-  REAL:: tau_ae2b_int(KLON,KLEV,nbands)   ! Intermediate computation of epaisseur optique aerosol
-  REAL:: piz_ae2b_int(KLON,KLEV,nbands)   ! Intermediate computation of Single scattering albedo
-  REAL:: cg_ae2b_int(KLON,KLEV,nbands)    ! Intermediate computation of Assymetry parameter
+  REAL:: DELTA(klon,klev), rh(klon,klev), H
+  REAL:: tau_ae2b_int   ! Intermediate computation of epaisseur optique aerosol
+  REAL:: piz_ae2b_int   ! Intermediate computation of Single scattering albedo
+  REAL:: cg_ae2b_int    ! Intermediate computation of Assymetry parameter
+  REAL :: Fact_RH(nbre_RH)
   REAL :: zrho
   REAL :: fac
@@ -146,4 +147,7 @@
   REAL:: piz_aeri_2bands(nbands,naero_insoluble)           !-- unit
 
+  INTEGER :: id
+  LOGICAL :: used_aer(naero_tot)
+  REAL :: tmp_var, tmp_var_pi
 
   DATA presnivs_19/&
@@ -586,25 +590,20 @@
 
   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
-        mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9
+    DO i=1, klon
+      zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
+!CDIR UNROLL=naero_spc
+      mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9
 !RAF zrho
-        mass_temp_pi(i,k,:) = m_allaer_pi(i,k,:) / zrho / 1.e+9
-     ENDDO
+!CDIR UNROLL=naero_spc
+      mass_temp_pi(i,k,:) = m_allaer_pi(i,k,:) / zrho / 1.e+9
+      zdp1(i,k)=pdel(i,k)/(gravit*delt)      ! air mass auxiliary  variable --> zdp1 [kg/(m^2 *s)]
+    ENDDO
   ENDDO
 
   IF (flag_aerosol .EQ. 1) THEN 
-     nb_aer = 1
+     nb_aer = 2
      ALLOCATE (aerosol_name(nb_aer)) 
      aerosol_name(1) = id_ASSO4M
-
+     aerosol_name(2) = id_CSSO4M
   ELSEIF (flag_aerosol .EQ. 2) THEN
      nb_aer = 2
@@ -650,23 +649,30 @@
   !      compute optical_thickness_at_gridpoint_per_species
 
-  tau_ae(:,:,:,:)=0.
-!RAF
-  tau_ae_pi(:,:,:,:)=0.
-  piz_ae(:,:,:,:)=0.
-  cg_ae(:,:,:,:)=0.
-  tau_allaer(:,:,:,:)=0.
-  piz_allaer(:,:,:,:)=0.
-  cg_allaer(:,:,:,:)=0.
-
-  !
-  ! Calculations that need to be done since we are not in the subroutines INCA
-  !      
-  ! air mass auxiliary  variable --> zdp1 [kg/(m^2 *s)]
-  zdp1(:,:)=pdel(:,:)/(gravit*delt) 
-
-
+
+
+!!CDIR ON_ADB(RH_tab)
+!CDIR ON_ADB(fact_RH)
+!CDIR SHORTLOOP
+  DO n=1,nbre_RH-1
+    fact_RH(n)=1./(RH_tab(n+1)-RH_tab(n))
+  ENDDO
+   
+  DO k=1, KLEV
+!!CDIR ON_ADB(RH_tab)
+!CDIR ON_ADB(fact_RH)
+    DO i=1, KLON
+      rh(i,k)=MIN(RHcl(i,k)*100.,RH_MAX)
+      RH_num(i,k) = INT( rh(i,k)/10. + 1.)
+      IF (rh(i,k).GT.85.) RH_num(i,k)=10
+      IF (rh(i,k).GT.90.) RH_num(i,k)=11
+      
+      DELTA(i,k)=(rh(i,k)-RH_tab(RH_num(i,k)))*fact_RH(RH_num(i,k))
+    ENDDO
+  ENDDO
+
+  used_aer(:)=.FALSE.
+    
   DO m=1,nb_aer   ! tau is only computed for each mass
-
-     fac=1.0
+    fac=1.0
      IF (aerosol_name(m).EQ.id_ASBCM) THEN
          soluble=.TRUE.
@@ -715,130 +721,243 @@
      ENDIF
 
-
-     tau_ae2b_int(:,:,:)=0.
-     piz_ae2b_int(:,:,:)=0.
-     cg_ae2b_int(:,:,:)=0.
-
-     DO inu=1,nbands
+    id=aerosol_name(m)
+    used_aer(id)=.TRUE.
+
+     
+    IF (soluble) THEN
+
+      IF (spss.NE.0) THEN
+
+         IF (spss.EQ.1) THEN !accumulation mode
+            DO k=1, KLEV
+!CDIR ON_ADB(A1_ASSSM_b1)
+!CDIR ON_ADB(A2_ASSSM_b1)
+!CDIR ON_ADB(A3_ASSSM_b1)
+!CDIR ON_ADB(B1_ASSSM_b1)
+!CDIR ON_ADB(B2_ASSSM_b1)
+!CDIR ON_ADB(C1_ASSSM_b1)
+!CDIR ON_ADB(C2_ASSSM_b2)
+!CDIR ON_ADB(A1_ASSSM_b2)
+!CDIR ON_ADB(A2_ASSSM_b2)
+!CDIR ON_ADB(A3_ASSSM_b2)
+!CDIR ON_ADB(B1_ASSSM_b2)
+!CDIR ON_ADB(B2_ASSSM_b2)
+!CDIR ON_ADB(C1_ASSSM_b2)
+!CDIR ON_ADB(C2_ASSSM_b2)
+              DO i=1, KLON
+                H=rh(i,k)/100
+                tmp_var=mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+                tmp_var_pi=mass_temp_pi(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+
+                ! band 1
+                tau_ae2b_int=A1_ASSSM_b1(k)+A2_ASSSM_b1(k)*H+A3_ASSSM_b1(k)/(H-1.05)
+                piz_ae2b_int=1-B1_ASSSM_b1(k)-B2_ASSSM_b1(k)*H
+                cg_ae2b_int=C1_ASSSM_b1(k)+C2_ASSSM_b1(k)*H
+
+                tau_ae(i,k,id,1) = tmp_var*tau_ae2b_int
+                tau_ae_pi(i,k,id,1) =  tmp_var_pi* tau_ae2b_int
+                piz_ae(i,k,id,1) = piz_ae2b_int
+                cg_ae(i,k,id,1)= cg_ae2b_int
+                
+                !band 2
+                tau_ae2b_int=A1_ASSSM_b2(k)+A2_ASSSM_b2(k)*H+A3_ASSSM_b2(k)/(H-1.05)
+                piz_ae2b_int=1-B1_ASSSM_b2(k)-B2_ASSSM_b2(k)*H
+                cg_ae2b_int=C1_ASSSM_b2(k)+C2_ASSSM_b2(k)*H
+
+                tau_ae(i,k,id,2) = tmp_var*tau_ae2b_int
+                tau_ae_pi(i,k,id,2) =  tmp_var_pi* tau_ae2b_int
+                piz_ae(i,k,id,2) = piz_ae2b_int
+                cg_ae(i,k,id,2)= cg_ae2b_int
+
+              ENDDO
+            ENDDO
+          ENDIF
+
+          IF (spss.EQ.2) THEN !coarse mode
+            DO k=1, KLEV
+!CDIR ON_ADB(A1_CSSSM_b1)
+!CDIR ON_ADB(A2_CSSSM_b1)
+!CDIR ON_ADB(A3_CSSSM_b1)
+!CDIR ON_ADB(B1_CSSSM_b1)
+!CDIR ON_ADB(B2_CSSSM_b1)
+!CDIR ON_ADB(C1_CSSSM_b1)
+!CDIR ON_ADB(C2_CSSSM_b2)
+!CDIR ON_ADB(A1_CSSSM_b2)
+!CDIR ON_ADB(A2_CSSSM_b2)
+!CDIR ON_ADB(A3_CSSSM_b2)
+!CDIR ON_ADB(B1_CSSSM_b2)
+!CDIR ON_ADB(B2_CSSSM_b2)
+!CDIR ON_ADB(C1_CSSSM_b2)
+!CDIR ON_ADB(C2_CSSSM_b2)
+              DO i=1, KLON
+                H=rh(i,k)/100
+                tmp_var=mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+                tmp_var_pi=mass_temp_pi(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+                ! band 1
+                tau_ae2b_int=A1_CSSSM_b1(k)+A2_CSSSM_b1(k)*H+A3_CSSSM_b1(k)/(H-1.05)
+                piz_ae2b_int=1-B1_CSSSM_b1(k)-B2_CSSSM_b1(k)*H
+                cg_ae2b_int=C1_CSSSM_b1(k)+C2_CSSSM_b1(k)*H
+
+                tau_ae(i,k,id,1) = tmp_var*tau_ae2b_int
+                tau_ae_pi(i,k,id,1) =  tmp_var_pi* tau_ae2b_int
+                piz_ae(i,k,id,1) = piz_ae2b_int
+                cg_ae(i,k,id,1)= cg_ae2b_int
+
+                ! band 2
+                tau_ae2b_int=A1_CSSSM_b2(k)+A2_CSSSM_b2(k)*H+A3_CSSSM_b2(k)/(H-1.05)
+                piz_ae2b_int=1-B1_CSSSM_b2(k)-B2_CSSSM_b2(k)*H
+                cg_ae2b_int=C1_CSSSM_b2(k)+C2_CSSSM_b2(k)*H
+
+                tau_ae(i,k,id,2) = tmp_var*tau_ae2b_int
+                tau_ae_pi(i,k,id,2) =  tmp_var_pi* tau_ae2b_int
+                piz_ae(i,k,id,2) = piz_ae2b_int
+                cg_ae(i,k,id,2)= cg_ae2b_int
+
+             ENDDO
+           ENDDO
+         ENDIF
+
+         IF (spss.EQ.3) THEN !super coarse mode
+            DO k=1, KLEV
+!CDIR ON_ADB(A1_SSSSM_b1)
+!CDIR ON_ADB(A2_SSSSM_b1)
+!CDIR ON_ADB(A3_SSSSM_b1)
+!CDIR ON_ADB(B1_SSSSM_b1)
+!CDIR ON_ADB(B2_SSSSM_b1)
+!CDIR ON_ADB(C1_SSSSM_b1)
+!CDIR ON_ADB(C2_SSSSM_b2)
+!CDIR ON_ADB(A1_SSSSM_b2)
+!CDIR ON_ADB(A2_SSSSM_b2)
+!CDIR ON_ADB(A3_SSSSM_b2)
+!CDIR ON_ADB(B1_SSSSM_b2)
+!CDIR ON_ADB(B2_SSSSM_b2)
+!CDIR ON_ADB(C1_SSSSM_b2)
+!CDIR ON_ADB(C2_SSSSM_b2)
+              DO i=1, KLON
+                H=rh(i,k)/100
+                tmp_var=mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+                tmp_var_pi=mass_temp_pi(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+
+                ! band 1 
+                tau_ae2b_int=A1_SSSSM_b1(k)+A2_SSSSM_b1(k)*H+A3_SSSSM_b1(k)/(H-1.05)
+                piz_ae2b_int=1-B1_SSSSM_b1(k)-B2_SSSSM_b1(k)*H
+                cg_ae2b_int=C1_SSSSM_b1(k)+C2_SSSSM_b1(k)*H
+
+                tau_ae(i,k,id,1) = tmp_var*tau_ae2b_int
+                tau_ae_pi(i,k,id,1) =  tmp_var_pi* tau_ae2b_int
+                piz_ae(i,k,id,1) = piz_ae2b_int
+                cg_ae(i,k,id,1)= cg_ae2b_int
+
+                ! band 2
+                tau_ae2b_int=A1_SSSSM_b2(k)+A2_SSSSM_b2(k)*H+A3_SSSSM_b2(k)/(H-1.05)
+                piz_ae2b_int=1-B1_SSSSM_b2(k)-B2_SSSSM_b2(k)*H
+                cg_ae2b_int=C1_SSSSM_b2(k)+C2_SSSSM_b2(k)*H
+
+                tau_ae(i,k,id,2) = tmp_var*tau_ae2b_int
+                tau_ae_pi(i,k,id,2) =  tmp_var_pi* tau_ae2b_int
+                piz_ae(i,k,id,2) = piz_ae2b_int
+                cg_ae(i,k,id,2)= cg_ae2b_int
+
+              ENDDO
+            ENDDO
+          ENDIF
+
+        ELSE
+                        
+!CDIR ON_ADB(alpha_aers_2bands)
+!CDIR ON_ADB(piz_aers_2bands)
+!CDIR ON_ADB(cg_aers_2bands)
+          DO k=1, KLEV
+            DO i=1, KLON
+              tmp_var=mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+              tmp_var_pi=mass_temp_pi(i,k,spsol)*1000.*zdp1(i,k)*delt*fac
+!CDIR UNROLL=nbands
+              DO inu=1,nbands
+
+                tau_ae2b_int= alpha_aers_2bands(RH_num(i,k),inu,spsol)+ & 
+                              DELTA(i,k)* (alpha_aers_2bands(RH_num(i,k)+1,inu,spsol) - & 
+                              alpha_aers_2bands(RH_num(i,k),inu,spsol))
+                      
+                piz_ae2b_int = piz_aers_2bands(RH_num(i,k),inu,spsol) + & 
+                               DELTA(i,k)* (piz_aers_2bands(RH_num(i,k)+1,inu,spsol) - & 
+                               piz_aers_2bands(RH_num(i,k),inu,spsol))
+                      
+                cg_ae2b_int = cg_aers_2bands(RH_num(i,k),inu,spsol) + & 
+                              DELTA(i,k)* (cg_aers_2bands(RH_num(i,k)+1,inu,spsol) - & 
+                              cg_aers_2bands(RH_num(i,k),inu,spsol))
+
+                tau_ae(i,k,id,inu) = tmp_var*tau_ae2b_int
+                tau_ae_pi(i,k,id,inu) =  tmp_var_pi* tau_ae2b_int
+                piz_ae(i,k,id,inu) = piz_ae2b_int
+                cg_ae(i,k,id,inu)= cg_ae2b_int
+                         
+              ENDDO
+            ENDDO
+          ENDDO
+        
+        ENDIF                     
+
+      ELSE                                                    ! For all aerosol insoluble components
+
+!CDIR ON_ADB(alpha_aers_2bands)
+!CDIR ON_ADB(piz_aers_2bands)
+!CDIR ON_ADB(cg_aers_2bands)
         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
-                      IF (spss.EQ.1) THEN !accumulation mode
-                          ! band 1
-                          tau_ae2b_int(i,k,1)=A1_ASSSM_b1(k)+A2_ASSSM_b1(k)*H+A3_ASSSM_b1(k)/(H-1.05)
-                          piz_ae2b_int(i,k,1)=1-B1_ASSSM_b1(k)-B2_ASSSM_b1(k)*H
-                          cg_ae2b_int(i,k,1)=C1_ASSSM_b1(k)+C2_ASSSM_b1(k)*H
-                          !band 2
-                          tau_ae2b_int(i,k,2)=A1_ASSSM_b2(k)+A2_ASSSM_b2(k)*H+A3_ASSSM_b2(k)/(H-1.05)
-                          piz_ae2b_int(i,k,2)=1-B1_ASSSM_b2(k)-B2_ASSSM_b2(k)*H
-                          cg_ae2b_int(i,k,2)=C1_ASSSM_b2(k)+C2_ASSSM_b2(k)*H
-                      ENDIF
-                      IF (spss.EQ.2) THEN !coarse mode
-                          ! band 1
-                          tau_ae2b_int(i,k,1)=A1_CSSSM_b1(k)+A2_CSSSM_b1(k)*H+A3_CSSSM_b1(k)/(H-1.05)
-                          piz_ae2b_int(i,k,1)=1-B1_CSSSM_b1(k)-B2_CSSSM_b1(k)*H
-                          cg_ae2b_int(i,k,1)=C1_CSSSM_b1(k)+C2_CSSSM_b1(k)*H
-                          ! band 2
-                          tau_ae2b_int(i,k,2)=A1_CSSSM_b2(k)+A2_CSSSM_b2(k)*H+A3_CSSSM_b2(k)/(H-1.05)
-                          piz_ae2b_int(i,k,2)=1-B1_CSSSM_b2(k)-B2_CSSSM_b2(k)*H
-                          cg_ae2b_int(i,k,2)=C1_CSSSM_b2(k)+C2_CSSSM_b2(k)*H
-                      ENDIF
-                      IF (spss.EQ.3) THEN !super coarse mode
-                          ! band 1 
-                          tau_ae2b_int(i,k,1)=A1_SSSSM_b1(k)+A2_SSSSM_b1(k)*H+A3_SSSSM_b1(k)/(H-1.05)
-                          piz_ae2b_int(i,k,1)=1-B1_SSSSM_b1(k)-B2_SSSSM_b1(k)*H
-                          cg_ae2b_int(i,k,1)=C1_SSSSM_b1(k)+C2_SSSSM_b1(k)*H
-                          ! band 2
-                          tau_ae2b_int(i,k,2)=A1_SSSSM_b2(k)+A2_SSSSM_b2(k)*H+A3_SSSSM_b2(k)/(H-1.05)
-                          piz_ae2b_int(i,k,2)=1-B1_SSSSM_b2(k)-B2_SSSSM_b2(k)*H
-                          cg_ae2b_int(i,k,2)=C1_SSSSM_b2(k)+C2_SSSSM_b2(k)*H
-                      ENDIF
-                  ELSE
-                      tau_ae2b_int(i,k,inu)= &
-                         alpha_aers_2bands(RH_num,inu,spsol)+ & 
-                         DELTA* (alpha_aers_2bands(RH_num+1,inu,spsol) - & 
-                         alpha_aers_2bands(RH_num,inu,spsol))
-                      
-                      piz_ae2b_int(i,k,inu)= &
-                         piz_aers_2bands(RH_num,inu,spsol) + & 
-                         DELTA* (piz_aers_2bands(RH_num+1,inu,spsol) - & 
-                         piz_aers_2bands(RH_num,inu,spsol))
-                      
-                      cg_ae2b_int(i,k,inu)= &
-                         cg_aers_2bands(RH_num,inu,spsol) + & 
-                         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
-!RAF fac
-                  tau_ae_pi(i,k,aerosol_name(m),inu) = &
-                     mass_temp_pi(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)
-                 cg_ae2b_int(i,k,inu) = cg_aeri_2bands(inu,spinsol) 
-
-                 tau_ae(i,k,aerosol_name(m),inu) = &
-                      mass_temp(i,k,naero_soluble+ spinsol)*1000.*zdp1(i,k)* &
-                                               delt*tau_ae2b_int(i,k,inu)*fac
-!RAF fac
-                 tau_ae_pi(i,k,aerosol_name(m),inu) = &
-                    mass_temp_pi(i,k,naero_soluble+ spinsol)*1000.*zdp1(i,k)* &
-                                                 delt*tau_ae2b_int(i,k,inu)*fac
-              ENDIF
-
-              piz_ae(i,k,aerosol_name(m),inu) = piz_ae2b_int(i,k,inu)
-
-              cg_ae(i,k,aerosol_name(m),inu)= cg_ae2b_int(i,k,inu)
-
-           ENDDO    ! nbands : boucle sur les bandes spectrale
-        ENDDO     ! klon : Boucle sur les points geographiques (grille horizontale)
-     ENDDO     ! klev : Boucle sur les niveaux verticaux
-  ENDDO     ! nb_aer : Boucle  sur les masses de traceurs
-
+          DO i=1, KLON
+            tmp_var=mass_temp(i,k,naero_soluble+ spinsol)*1000.*zdp1(i,k)*delt*fac
+            tmp_var_pi=mass_temp_pi(i,k,naero_soluble+spinsol)*1000.*zdp1(i,k)*delt*fac
+!CDIR UNROLL=nbands
+            DO inu=1,nbands
+              tau_ae2b_int = alpha_aeri_2bands(inu,spinsol)
+              piz_ae2b_int = piz_aeri_2bands(inu,spinsol)
+              cg_ae2b_int = cg_aeri_2bands(inu,spinsol) 
+
+              tau_ae(i,k,id,inu) = tmp_var*tau_ae2b_int
+              tau_ae_pi(i,k,id,inu) = tmp_var_pi*tau_ae2b_int
+              piz_ae(i,k,id,inu) = piz_ae2b_int
+              cg_ae(i,k,id,inu)= cg_ae2b_int
+            ENDDO
+          ENDDO
+        ENDDO
+
+      ENDIF ! soluble
+
+    ENDDO  ! nb_aer  
+
+  DO m=1,nb_aer   
+    IF (.NOT. used_aer(m)) THEN
+      tau_ae(:,:,:,:)=0.
+      tau_ae_pi(:,:,:,:)=0.
+      piz_ae(:,:,:,:)=0.
+      cg_ae(:,:,:,:)=0.
+    ENDIF
+  ENDDO
 
   DO inu=1, nbands
-     DO mrfspecies=1,naero_grp
+    DO mrfspecies=1,naero_grp
+      IF (mrfspecies .EQ. 2) THEN             ! = total aerosol AER	 
         DO k=1, KLEV
-           DO i=1, KLON
-              IF (mrfspecies .EQ. 2) THEN             ! = total aerosol AER	 
-	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu)+ &
-                      tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu)+ &						     
-                      tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu)+ &	
-                      tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+tau_ae(i,k,id_SSSSM,inu)+ &
-                      tau_ae(i,k,id_CIDUSTM,inu)
-	         tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
+          DO i=1, KLON
+            tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu)+ &
+                                           tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu)+   &						     
+                                           tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu)+ &	
+                                           tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+   &
+                                           tau_ae(i,k,id_SSSSM,inu)+ tau_ae(i,k,id_CIDUSTM,inu)
+	     tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
                  
-	 	 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)+ &
-                      tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)+ &
-                      tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)+ &
-                      tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)+ &
-                      tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)+ &
-                      tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)+ &	
-                      tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)+ &
-                      tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)+ &
-                      tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)+ &
-                      tau_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu))/tau_allaer(i,k,mrfspecies,inu)
-	         piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
-
-	 	 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ &
+             piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)+ &
+                                             tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)+ &
+                                             tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)+ &
+                                             tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)+ &
+                                             tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)+ &
+                                             tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)+ &	
+                                             tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)+ &
+                                             tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)+ &
+                                             tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)+ &
+                                             tau_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu)) &
+                                            /tau_allaer(i,k,mrfspecies,inu)
+	     piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
+
+             cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ &
                       tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)*cg_ae(i,k,id_CSSO4M,inu)+ &
                       tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)*cg_ae(i,k,id_ASBCM,inu)+ &
@@ -851,6 +970,11 @@
                       tau_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu)*cg_ae(i,k,id_CIDUSTM,inu))/ &
                       (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
-
-              ELSEIF (mrfspecies .EQ. 3) THEN             ! = natural aerosol NAT
+          ENDDO    
+        ENDDO 
+
+      ELSEIF (mrfspecies .EQ. 3) THEN             ! = natural aerosol NAT
+
+        DO k=1, KLEV
+          DO i=1, KLON
 !RAF
 	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae_pi(i,k,id_ASSO4M,inu)+ &
@@ -892,66 +1016,106 @@
                       cg_ae(i,k,id_CIDUSTM,inu))/ &
                       (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
-
-!END RAF
-              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)
-	         tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
-		 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) &
+          ENDDO
+        ENDDO
+                   
+      ELSEIF (mrfspecies .EQ. 4) THEN             ! = BC
+        DO k=1, KLEV
+          DO i=1, KLON
+	    tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu)
+	    tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
+	    piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) &
                       +tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu))/ &
                       tau_allaer(i,k,mrfspecies,inu)
-	         piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
-		 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) *cg_ae(i,k,id_ASBCM,inu)&
+	    piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
+            cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) *cg_ae(i,k,id_ASBCM,inu)&
                       +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)
-	         tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
-		 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) &
+          ENDDO
+        ENDDO
+              
+      ELSEIF (mrfspecies .EQ. 5) THEN             ! = SO4
+
+        DO k=1, KLEV
+          DO i=1, KLON
+            tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu)
+	    tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
+            piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) &
                       +tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu))/ &
                       tau_allaer(i,k,mrfspecies,inu)
-	         piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
-		 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) *cg_ae(i,k,id_CSSO4M,inu)&
+	    piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
+            cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) *cg_ae(i,k,id_CSSO4M,inu)&
                       +tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu))/ &
                       (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
-
-              ELSEIF (mrfspecies .EQ. 6) THEN             ! = POM
-	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu)
-	         tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
-		 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) &
+          ENDDO
+        ENDDO
+
+      ELSEIF (mrfspecies .EQ. 6) THEN             ! = POM
+
+        DO k=1, KLEV
+          DO i=1, KLON
+            tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu)
+            tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
+	    piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) &
                       +tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu))/ &
                       tau_allaer(i,k,mrfspecies,inu)
-	         piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
-		 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) *cg_ae(i,k,id_ASPOMM,inu)&
+	    piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
+	    cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) *cg_ae(i,k,id_ASPOMM,inu)&
                       +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)
-	         tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
-		 piz_allaer(i,k,mrfspecies,inu)=piz_ae(i,k,id_CIDUSTM,inu)
-		 cg_allaer(i,k,mrfspecies,inu)=cg_ae(i,k,id_CIDUSTM,inu)
-
-              ELSEIF (mrfspecies .EQ. 8) THEN             ! = SS
-	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+tau_ae(i,k,id_SSSSM,inu)
-	         tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
-		 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) &
-                      +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu) &
-                      +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu))/ &
-                      tau_allaer(i,k,mrfspecies,inu)
-	         piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
-		 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) *cg_ae(i,k,id_ASSSM,inu)&
-                      +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu) &
-                      +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu))/ &
-                      (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
-
-              ELSEIF (mrfspecies .EQ. 9) THEN             ! = NO3
-	 	 tau_allaer(i,k,mrfspecies,inu)=0.   ! preliminary
-		 piz_allaer(i,k,mrfspecies,inu)=0.
-		 cg_allaer(i,k,mrfspecies,inu)=0.
-              ENDIF
-           ENDDO
+          ENDDO
         ENDDO
-     ENDDO
+              
+      ELSEIF (mrfspecies .EQ. 7) THEN             ! = DUST
+
+        DO k=1, KLEV
+          DO i=1, KLON
+            tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_CIDUSTM,inu)
+	    tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
+            piz_allaer(i,k,mrfspecies,inu)=piz_ae(i,k,id_CIDUSTM,inu)
+	    cg_allaer(i,k,mrfspecies,inu)=cg_ae(i,k,id_CIDUSTM,inu)
+          ENDDO
+        ENDDO
+
+      ELSEIF (mrfspecies .EQ. 8) THEN             ! = SS
+
+        DO k=1, KLEV
+          DO i=1, KLON
+            tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+tau_ae(i,k,id_SSSSM,inu)
+            tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20)
+            piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) &
+                    +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu) &
+                    +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu))/ &
+                    tau_allaer(i,k,mrfspecies,inu)
+            piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20)
+            cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) *cg_ae(i,k,id_ASSSM,inu)&
+                    +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu) &
+                    +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu))/ &
+                    (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
+          ENDDO
+        ENDDO
+      
+      ELSEIF (mrfspecies .EQ. 9) THEN             ! = NO3
+      
+        DO k=1, KLEV
+          DO i=1, KLON
+            tau_allaer(i,k,mrfspecies,inu)=0.   ! preliminary
+            piz_allaer(i,k,mrfspecies,inu)=0.
+            cg_allaer(i,k,mrfspecies,inu)=0.
+          ENDDO
+        ENDDO
+      
+      ELSE
+
+        DO k=1, KLEV
+          DO i=1, KLON
+            tau_allaer(i,k,mrfspecies,inu)=0.  
+            piz_allaer(i,k,mrfspecies,inu)=0.
+            cg_allaer(i,k,mrfspecies,inu)=0.
+          ENDDO
+        ENDDO
+           
+      ENDIF
+
+    ENDDO
   ENDDO
 
