Index: LMDZ5/trunk/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 2004)
+++ LMDZ5/trunk/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 2004)
@@ -0,0 +1,413 @@
+!
+! $Id: aeropt_5wv.F90 1716 2013-01-22 13:36:30Z acozic $
+!
+
+SUBROUTINE AEROPT_5WV_RRTM(&
+   pdel, m_allaer, delt, &
+   RHcl, ai, flag_aerosol, &
+   pplay, t_seri, &
+   tausum, tau )
+
+  USE DIMPHY
+  USE aero_mod
+  USE phys_local_var_mod, only: od550aer,od865aer,ec550aer,od550lt1aer
+
+  !
+  !    Yves Balkanski le 12 avril 2006
+  !    Celine Deandreis
+  !    Anne Cozic  Avril 2009
+  !    a partir d'une sous-routine de Johannes Quaas pour les sulfates
+  !    Olivier Boucher mars 2014 pour adaptation RRTM
+  !
+  ! Refractive indices for seasalt come from Shettle and Fenn (1979)
+  !
+  ! Refractive indices from water come from Hale and Querry (1973)
+  !
+  ! Refractive indices from Ammonium Sulfate Toon and Pollack (1976)
+  !
+  ! Refractive indices for Dust, internal mixture of minerals coated with 1.5% hematite 
+  ! by Volume (Balkanski et al., 2006)
+  !
+  ! Refractive indices for POM: Kinne (pers. Communication 
+  !
+  ! Refractive index for BC from Shettle and Fenn (1979)
+  !
+  ! Shettle, E. P., & Fenn, R. W. (1979), Models for the aerosols of the lower atmosphere and 
+  ! the effects of humidity variations on their optical properties, U.S. Air Force Geophysics 
+  ! Laboratory Rept. AFGL-TR-79-0214, Hanscomb Air Force Base, MA.
+  !
+  ! Hale, G. M. and M. R. Querry, Optical constants of water in the 200-nm to 200-m 
+  ! wavelength region, Appl. Opt., 12, 555-563, 1973.
+  !
+  ! Toon, O. B. and J. B. Pollack, The optical constants of several atmospheric aerosol species:
+  ! Ammonium sulfate, aluminum oxide, and sodium chloride, J. Geohys. Res., 81, 5733-5748,
+  ! 1976.
+  !
+  ! Balkanski, Y., M. Schulz, T. Claquin And O. Boucher, Reevaluation of mineral aerosol 
+  ! radiative forcings suggests a better agreement with satellite and AERONET data, Atmospheric 
+  ! Chemistry and Physics Discussions., 6, pp 8383-8419, 2006.
+  !
+  IMPLICIT NONE
+  INCLUDE "YOMCST.h"
+  !
+  ! Input arguments:
+  !
+  REAL, DIMENSION(klon,klev), INTENT(in)   :: pdel
+  REAL, INTENT(in)                         :: delt
+  REAL, DIMENSION(klon,klev,naero_spc), INTENT(in) :: m_allaer
+  REAL, DIMENSION(klon,klev), INTENT(in)   :: RHcl     ! humidite relative ciel clair
+  INTEGER,INTENT(in)                       :: flag_aerosol
+  REAL, DIMENSION(klon,klev), INTENT(in)   :: pplay
+  REAL, DIMENSION(klon,klev), INTENT(in)   :: t_seri
+  !
+  ! Output arguments:
+  !
+  REAL, DIMENSION(klon), INTENT(out)                      :: ai      ! POLDER aerosol index 
+  REAL, DIMENSION(klon,nwave,naero_tot), INTENT(out)      :: tausum
+  REAL, DIMENSION(klon,klev,nwave,naero_tot), INTENT(out) :: tau
+
+  !
+  ! Local
+  !
+  INTEGER, PARAMETER :: las = nwave
+  LOGICAL :: soluble
+  
+  INTEGER :: i, k, m
+  INTEGER :: spsol, spinsol, la
+  INTEGER :: RH_num(klon,klev)
+  INTEGER, PARAMETER :: la443 = 1
+  INTEGER, PARAMETER :: la550 = 2
+  INTEGER, PARAMETER :: la670 = 3
+  INTEGER, PARAMETER :: la765 = 4
+  INTEGER, PARAMETER :: la865 = 5
+  INTEGER, PARAMETER :: nbre_RH=12
+  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
+
+  REAL :: zrho
+
+  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(klon,klev), rh(klon,klev)
+  REAL :: tau_ae5wv_int   ! Intermediate computation of epaisseur optique aerosol
+  REAL :: od670aer(klon)  ! epaisseur optique aerosol extinction 670 nm
+  REAL :: fac
+  REAL :: zdp1(klon,klev) 
+  INTEGER, ALLOCATABLE, DIMENSION(:)  :: aerosol_name
+  INTEGER :: nb_aer
+  
+  REAL :: dh(KLON,KLEV)
+  
+   ! Soluble components 1- BC soluble; 2- POM soluble; 3- SO4 acc.; 4- SO4 coarse; 5 seasalt super-coarse; 6 seasalt coarse; 7 seasalt acc.
+  REAL :: alpha_aers_5wv(nbre_RH,las,naero_soluble)   ! Ext. coeff. ** m2/g 
+   ! Insoluble components 1- Dust: 2- BC; 3- POM
+  REAL :: alpha_aeri_5wv(las,naero_insoluble)         ! Ext. coeff. ** m2/g 
+
+  REAL, DIMENSION(klon,klev,naero_spc) :: mass_temp
+  
+  !
+  ! Proprietes optiques
+  !
+  REAL :: fact_RH(nbre_RH)
+  LOGICAL :: used_tau(naero_tot)
+  INTEGER :: n
+  
+! From here on we look at the optical parameters at 5 wavelengths:  
+! 443nm, 550, 670, 765 and 865 nm 
+!                                   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 (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 seasalt Super Coarse Soluble (SS)
+        0.218, 0.240, 0.262, 0.311, 0.391, 0.452, &
+        0.544, 0.614, 0.886, 1.018, 1.260, 2.336, &
+        0.221, 0.243, 0.266, 0.315, 0.396, 0.455, &
+        0.551, 0.622, 0.896, 1.029, 1.270, 2.352, &
+        0.224, 0.246, 0.270, 0.319, 0.400, 0.461, &
+        0.558, 0.626, 0.904, 1.039, 1.278, 2.362, &
+        0.227, 0.249, 0.272, 0.322, 0.404, 0.466, &
+        0.560, 0.632, 0.907, 1.041, 1.284, 2.375, &
+        0.230, 0.252, 0.276, 0.325, 0.408, 0.469, &
+        0.567, 0.636, 0.917, 1.052, 1.293, 2.383, &
+                           ! seasalt seasalt Coarse Soluble (CS)      
+        0.578, 0.632, 0.691, 0.814, 1.017, 1.171, &
+        1.407, 1.585, 2.264, 2.594, 3.197, 5.882, &
+        0.598, 0.654, 0.711, 0.836, 1.041, 1.201, &
+        1.445, 1.622, 2.313, 2.647, 3.256, 5.961, &
+        0.619, 0.676, 0.736, 0.862, 1.072, 1.234, &
+        1.484, 1.659, 2.361, 2.698, 3.304, 6.014, &
+        0.633, 0.692, 0.754, 0.884, 1.095, 1.262, &
+        1.510, 1.693, 2.395, 2.734, 3.352, 6.090, &
+        0.648, 0.708, 0.771, 0.902, 1.119, 1.288, &
+        1.545, 1.725, 2.441, 2.784, 3.402, 6.144, &
+                           ! seasalt seasalt Accumulation Soluble (AS)
+        4.432, 5.160, 5.940, 7.646,10.520,11.464, &
+       14.686,16.949,25.977,30.118,37.327,65.333, &
+        3.298, 3.916, 4.594, 6.127, 8.845, 9.446, &
+       12.511,14.710,24.073,28.546,36.578,69.050, &
+        2.340, 2.826, 3.371, 4.638, 6.992, 7.360, &
+       10.047,12.016,20.909,25.349,33.605,69.554, &
+        1.789, 2.185, 2.633, 3.698, 5.733, 5.994, &
+        8.345,10.094,18.293,22.511,30.548,67.717, &
+        1.359, 1.676, 2.040, 2.916, 4.637, 4.829, &
+        6.843, 8.364,15.723,19.614,27.195,64.359 /
+
+  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/ 
+  ! 
+  ! Initialisations
+  ai(:) = 0.
+  tausum(:,:,:) = 0.
+
+  DO k=1, klev
+    DO i=1, klon
+      zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
+      dh(i,k)=pdel(i,k)/(RG*zrho)
+!CDIR UNROLL=naero_spc
+      mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9
+      zdp1(i,k)=pdel(i,k)/(RG*delt)     ! air mass auxiliary  variable --> zdp1 [kg/(m^2 *s)]
+    ENDDO
+  ENDDO
+
+  IF (flag_aerosol .EQ. 1) THEN 
+     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
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASBCM
+     aerosol_name(2) = id_AIBCM
+  ELSEIF (flag_aerosol .EQ. 3) THEN 
+     nb_aer = 2
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASPOMM
+     aerosol_name(2) = id_AIPOMM
+  ELSEIF (flag_aerosol .EQ. 4) THEN 
+     nb_aer = 3
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_CSSSM
+     aerosol_name(2) = id_SSSSM
+     aerosol_name(3) = id_ASSSM
+  ELSEIF (flag_aerosol .EQ. 5) THEN 
+     nb_aer = 1
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_CIDUSTM
+  ELSEIF (flag_aerosol .EQ. 6) THEN 
+     nb_aer = 10
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASSO4M      
+     aerosol_name(2) = id_ASBCM
+     aerosol_name(3) = id_AIBCM
+     aerosol_name(4) = id_ASPOMM
+     aerosol_name(5) = id_AIPOMM
+     aerosol_name(6) = id_CSSSM
+     aerosol_name(7) = id_SSSSM
+     aerosol_name(8) = id_ASSSM
+     aerosol_name(9) = id_CIDUSTM
+     aerosol_name(10) = id_CSSO4M
+  ENDIF
+
+  ! 
+  ! Loop over modes, use of precalculated nmd and corresponding sigma
+  !    loop over wavelengths
+  !    for each mass species in mode
+  !      interpolate from Sext to retrieve Sext_at_gridpoint_per_species
+  !      compute optical_thickness_at_gridpoint_per_species
+  !
+  ! Calculations that need to be done since we are not in the subroutines INCA
+  !      
+
+!CDIR ON_ADB(RH_tab)
+!CDIR ON_ADB(fact_RH)
+!CDIR NOVECTOR
+  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
+
+!CDIR SHORTLOOP  
+  used_tau(:)=.FALSE.
+    
+  DO m=1,nb_aer   ! tau is only computed for each mass    
+    fac=1.0
+    IF (aerosol_name(m).EQ.id_ASBCM) THEN
+        soluble=.TRUE.
+        spsol=1
+    ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN 
+        soluble=.TRUE.
+        spsol=2 
+    ELSEIF (aerosol_name(m).EQ.id_ASSO4M) THEN
+        soluble=.TRUE.
+        spsol=3
+        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
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+    ELSEIF (aerosol_name(m).EQ.id_SSSSM) THEN 
+        soluble=.TRUE.
+        spsol=5
+    ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN 
+        soluble=.TRUE.
+        spsol=6
+    ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN
+        soluble=.TRUE.
+        spsol=7
+    ELSEIF (aerosol_name(m).EQ.id_CIDUSTM) THEN 
+        soluble=.FALSE.
+        spinsol=1
+    ELSEIF  (aerosol_name(m).EQ.id_AIBCM) THEN 
+        soluble=.FALSE.
+        spinsol=2
+    ELSEIF (aerosol_name(m).EQ.id_AIPOMM) THEN 
+        soluble=.FALSE.
+        spinsol=3
+    ELSE 
+        CYCLE
+    ENDIF
+
+    IF (soluble) then
+      used_tau(spsol)=.TRUE.
+    ELSE
+      used_tau(naero_soluble+spinsol)=.TRUE.
+    ENDIF
+
+    DO la=1,las
+
+      IF (soluble) THEN            ! For soluble aerosol
+
+          DO k=1, KLEV
+            DO i=1, KLON
+              tau_ae5wv_int = alpha_aers_5wv(RH_num(i,k),la,spsol)+DELTA(i,k)* &
+                             (alpha_aers_5wv(RH_num(i,k)+1,la,spsol) - & 
+                              alpha_aers_5wv(RH_num(i,k),la,spsol))
+              tau(i,k,la,spsol) = mass_temp(i,k,spsol)*1000.*zdp1(i,k)   &
+                                 *tau_ae5wv_int*delt*fac
+              tausum(i,la,spsol)=tausum(i,la,spsol)+tau(i,k,la,spsol)
+            ENDDO
+          ENDDO
+
+      ELSE                         ! For insoluble aerosol
+
+        DO k=1, KLEV
+          DO i=1, KLON
+            tau_ae5wv_int = alpha_aeri_5wv(la,spinsol)
+            tau(i,k,la,naero_soluble+spinsol) = mass_temp(i,k,naero_soluble+spinsol)*1000.*zdp1(i,k)* &
+                                                tau_ae5wv_int*delt*fac
+            tausum(i,la,naero_soluble+spinsol)= tausum(i,la,naero_soluble+spinsol)  &
+                                               +tau(i,k,la,naero_soluble+spinsol)
+          ENDDO
+        ENDDO
+
+      ENDIF
+
+    ENDDO   ! Boucle sur les longueurs d'onde
+  ENDDO     ! Boucle sur les masses de traceurs
+
+  DO m=1,naero_tot
+    IF (.NOT.used_tau(m)) tau(:,:,:,m)=0.
+  ENDDO  
+
+  DO i=1, klon
+     od550aer(i)=0.
+     DO m=1,naero_tot
+        od550aer(i)=od550aer(i)+tausum(i,la550,m)
+     END DO
+  END DO
+
+  DO i=1, klon
+     od670aer(i)=0.
+     DO m=1,naero_tot
+        od670aer(i)=od670aer(i)+tausum(i,la670,m)
+     END DO
+  END DO
+
+  DO i=1, klon
+     od865aer(i)=0.
+     DO m=1,naero_tot
+        od865aer(i)=od865aer(i)+tausum(i,la865,m)
+     END DO
+  END DO
+
+  DO i=1, klon
+     DO k=1, KLEV
+        ec550aer(i,k)=0.
+        DO m=1,naero_tot
+           ec550aer(i,k)=ec550aer(i,k)+tau(i,k,la550,m)/dh(i,k)
+        END DO
+     END DO
+  END DO
+  
+  DO i=1, klon
+    ai(i)=-LOG(MAX(od670aer(i),1.e-8)/MAX(od865aer(i),1.e-8))/LOG(670./865.)
+  ENDDO
+
+  od550lt1aer(:)=tausum(:,la550,id_ASSO4M)+tausum(:,la550,id_ASBCM) +tausum(:,la550,id_AIBCM)+ &
+                 tausum(:,la550,id_ASPOMM)+tausum(:,la550,id_AIPOMM)+tausum(:,la550,id_ASSSM)+ &
+                 0.03*tausum(:,la550,id_CSSSM)+0.4*tausum(:,la550,id_CIDUSTM)
+
+  DEALLOCATE(aerosol_name) 
+  
+END SUBROUTINE AEROPT_5WV_RRTM
Index: LMDZ5/trunk/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 2004)
+++ LMDZ5/trunk/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 2004)
@@ -0,0 +1,550 @@
+!
+! $Id: aeropt_6bands.F90 1716 2013-01-22 13:36:30Z acozic $
+!
+SUBROUTINE AEROPT_6BANDS_RRTM ( &
+     pdel, m_allaer, delt, RHcl, &
+     tau_allaer, piz_allaer, &
+     cg_allaer, m_allaer_pi, &
+     flag_aerosol, pplay, t_seri )
+
+  USE dimphy
+  USE aero_mod
+  USE phys_local_var_mod, only: absvisaer
+
+  !    Yves Balkanski le 12 avril 2006
+  !    Celine Deandreis
+  !    Anne Cozic Avril 2009
+  !    a partir d'une sous-routine de Johannes Quaas pour les sulfates
+  !    Olivier Boucher février 2014 pour passage à RRTM
+  !    a partir des propriétés optiques fournies par Yves Balkanski
+  !
+  IMPLICIT NONE
+
+  INCLUDE "YOMCST.h"
+  INCLUDE "iniprint.h"
+  INCLUDE "clesphys.h"
+
+  !
+  ! Input arguments:
+  !
+  REAL, DIMENSION(klon,klev),     INTENT(in)  :: pdel
+  REAL,                           INTENT(in)  :: delt
+  REAL, DIMENSION(klon,klev,naero_spc),   INTENT(in)  :: m_allaer
+  REAL, DIMENSION(klon,klev,naero_spc),   INTENT(in)  :: m_allaer_pi
+  REAL, DIMENSION(klon,klev),     INTENT(in)  :: RHcl       ! humidite relative ciel clair
+  INTEGER,                        INTENT(in)  :: flag_aerosol
+  REAL, DIMENSION(klon,klev),     INTENT(in)  :: pplay
+  REAL, DIMENSION(klon,klev),     INTENT(in)  :: t_seri
+  !
+  ! Output arguments:
+  ! 1= total aerosols 
+  ! 2= natural aerosols
+  !
+  REAL, DIMENSION(klon,klev,2,nbands_rrtm), INTENT(out) :: tau_allaer ! epaisseur optique aerosol
+  REAL, DIMENSION(klon,klev,2,nbands_rrtm), INTENT(out) :: piz_allaer ! single scattering albedo aerosol
+  REAL, DIMENSION(klon,klev,2,nbands_rrtm), INTENT(out) :: cg_allaer  ! asymmetry parameter aerosol
+
+  !
+  ! Local
+  !
+  LOGICAL ::  soluble
+  INTEGER :: i, k,n, inu, m
+  INTEGER :: spsol, spinsol
+  INTEGER :: RH_num(klon,klev)
+
+  INTEGER, PARAMETER :: nb_level=19 ! number of vertical levels in DATA
+
+  INTEGER, PARAMETER :: naero_soluble=7    ! 1- BC soluble; 2- POM soluble; 3- SO4. acc. 4- SO4 coarse
+                                           ! 5- seasalt super coarse  6- seasalt coarse   7- seasalt acc.
+  INTEGER, PARAMETER :: naero_insoluble=3  ! 1- Dust; 2- BC insoluble; 3- POM insoluble
+
+  INTEGER, PARAMETER :: nbre_RH=12
+  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(klon,klev), rh(klon,klev)
+  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
+  REAL :: zdp1(klon,klev) 
+  INTEGER, ALLOCATABLE, DIMENSION(:)   :: aerosol_name
+  INTEGER :: nb_aer
+
+  REAL, DIMENSION(klon,klev,naero_spc) :: mass_temp
+  REAL, DIMENSION(klon,klev,naero_spc) :: mass_temp_pi
+  REAL, DIMENSION(klon,klev,naero_tot,nbands_rrtm) ::  tau_ae
+  REAL, DIMENSION(klon,klev,naero_tot,nbands_rrtm) ::  tau_ae_pi
+  REAL, DIMENSION(klon,klev,naero_tot,nbands_rrtm) ::  piz_ae
+  REAL, DIMENSION(klon,klev,naero_tot,nbands_rrtm) ::  cg_ae
+
+
+  !
+  ! Proprietes optiques
+  !
+  REAL:: alpha_aers_6bands(nbre_RH,nbands_rrtm,naero_soluble)   !--unit m2/g SO4
+  REAL:: alpha_aeri_6bands(nbands_rrtm,naero_insoluble)
+  REAL:: cg_aers_6bands(nbre_RH,nbands_rrtm,naero_soluble)      !--unit 
+  REAL:: cg_aeri_6bands(nbands_rrtm,naero_insoluble)
+  REAL:: piz_aers_6bands(nbre_RH,nbands_rrtm,naero_soluble)     !-- unit
+  REAL:: piz_aeri_6bands(nbands_rrtm,naero_insoluble)        !-- unit
+
+  INTEGER :: id
+  LOGICAL :: used_aer(naero_tot)
+  REAL :: tmp_var, tmp_var_pi
+
+!***************************************************************************
+!--the order of the soluble   species has to follow the spsol   index below
+!--the order of the insoluble species has to follow the spinsol index below
+
+  DATA alpha_aers_6bands/  & 
+       ! bc soluble AS
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
+  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
+  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
+  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
+  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
+       ! pom soluble AS
+  6.443, 6.443, 6.443, 6.443, 6.443, 7.100, 7.809, 9.277,10.721,10.882,11.056,12.628, &
+  6.443, 6.443, 6.443, 6.443, 6.443, 7.100, 7.809, 9.277,10.721,10.882,11.056,12.628, &
+  4.381, 4.381, 4.381, 4.381, 4.381, 4.828, 5.310, 6.309, 7.290, 7.400, 7.518, 8.587, &
+  1.846, 1.846, 1.846, 1.846, 1.846, 2.034, 2.237, 2.658, 3.072, 3.118, 3.168, 3.618, &
+  0.377, 0.377, 0.377, 0.377, 0.377, 0.415, 0.456, 0.542, 0.627, 0.636, 0.646, 0.738, &
+  0.052, 0.052, 0.052, 0.052, 0.052, 0.057, 0.063, 0.075, 0.087, 0.088, 0.089, 0.102, &
+       ! sulfate AS   
+  6.554, 6.554, 6.554, 7.223, 7.931, 8.665, 9.438,10.736,14.275,17.755,17.755,31.722, &
+  6.554, 6.554, 6.554, 7.223, 7.931, 8.665, 9.438,10.736,14.275,17.755,17.755,31.722, &
+  4.381, 4.381, 4.381, 4.828, 5.301, 5.792, 6.309, 7.176, 9.542,11.868,11.868,21.204, &
+  1.727, 1.727, 1.727, 1.903, 2.090, 2.283, 2.487, 2.829, 3.762, 4.679, 4.679, 8.359, &
+  0.312, 0.312, 0.312, 0.344, 0.378, 0.413, 0.450, 0.511, 0.680, 0.846, 0.846, 1.511, &
+  0.121, 0.121, 0.121, 0.134, 0.147, 0.161, 0.175, 0.199, 0.264, 0.329, 0.329, 0.588, &
+       ! sulfate coarse CS
+  0.693, 0.693, 0.693, 0.764, 0.839, 0.917, 0.999, 1.136, 1.510, 1.879, 1.879, 3.356, &
+  0.693, 0.693, 0.693, 0.764, 0.839, 0.917, 0.999, 1.136, 1.510, 1.879, 1.879, 3.356, &
+  0.715, 0.715, 0.715, 0.788, 0.865, 0.945, 1.029, 1.171, 1.557, 1.936, 1.936, 3.459, &
+  0.736, 0.736, 0.736, 0.811, 0.891, 0.973, 1.060, 1.206, 1.603, 1.994, 1.994, 3.563, &
+  0.711, 0.711, 0.711, 0.783, 0.860, 0.939, 1.023, 1.164, 1.548, 1.925, 1.925, 3.439, &
+  0.602, 0.602, 0.602, 0.664, 0.729, 0.796, 0.867, 0.986, 1.312, 1.631, 1.631, 2.915, &
+       ! seasalt seasalt Super Coarse Soluble (SS)
+  0.214, 0.235, 0.257, 0.305, 0.385, 0.444, 0.537, 0.606, 0.876, 1.006, 1.243, 2.313, &
+  0.217, 0.238, 0.261, 0.309, 0.389, 0.448, 0.543, 0.612, 0.882, 1.015, 1.255, 2.327, &
+  0.221, 0.243, 0.266, 0.315, 0.396, 0.456, 0.552, 0.621, 0.895, 1.028, 1.269, 2.351, &
+  0.230, 0.253, 0.276, 0.326, 0.408, 0.469, 0.568, 0.638, 0.918, 1.053, 1.298, 2.389, &
+  0.250, 0.273, 0.298, 0.349, 0.435, 0.501, 0.601, 0.674, 0.960, 1.098, 1.349, 2.466, &
+  0.279, 0.307, 0.336, 0.396, 0.493, 0.566, 0.679, 0.760, 1.072, 1.221, 1.486, 2.652, &
+       ! seasalt seasalt Coarse Soluble (CS)      
+  0.550, 0.605, 0.661, 0.782, 0.980, 1.131, 1.362, 1.536, 2.204, 2.528, 3.115, 5.773, &
+  0.568, 0.622, 0.681, 0.802, 1.004, 1.155, 1.393, 1.568, 2.246, 2.574, 3.174, 5.846, &
+  0.599, 0.655, 0.714, 0.838, 1.044, 1.202, 1.446, 1.621, 2.312, 2.645, 3.249, 5.950, &
+  0.647, 0.708, 0.771, 0.904, 1.121, 1.290, 1.546, 1.731, 2.447, 2.791, 3.416, 6.166, &
+  0.663, 0.735, 0.809, 0.964, 1.215, 1.386, 1.677, 1.883, 2.681, 3.059, 3.733, 6.647, &
+  0.479, 0.549, 0.623, 0.787, 1.065, 1.174, 1.491, 1.717, 2.648, 3.094, 3.901, 7.382, &
+       ! seasalt seasalt Accumulation Soluble (AS)
+  6.128, 6.867, 7.624, 9.184,11.602,13.194,15.955,17.858,24.856,28.034,33.562,56.971, &
+  5.080, 5.837, 6.636, 8.347,11.144,12.321,15.463,17.649,26.080,29.896,36.491,62.414, &
+  3.290, 3.899, 4.567, 6.075, 8.744, 9.355,12.369,14.532,23.747,28.159,36.101,68.531, &
+  1.389, 1.706, 2.070, 2.942, 4.643, 4.846, 6.836, 8.334,15.546,19.346,26.738,62.978, &
+  0.309, 0.394, 0.495, 0.754, 1.312, 1.377, 2.087, 2.654, 5.732, 7.538,11.368,35.212, &
+  0.021, 0.028, 0.036, 0.059, 0.114, 0.126, 0.208, 0.279, 0.718, 1.009, 1.691, 7.271  /
+
+  DATA alpha_aeri_6bands/  & 
+       ! dust insoluble CI
+  0.751, 0.751, 0.769, 0.772, 0.672, 0.437, &
+       ! bc insoluble AI
+  6.497, 6.497, 5.900, 4.284, 2.163, 0.966, &
+       ! pom insoluble AI
+  6.443, 6.443, 4.381, 1.846, 0.377, 0.052 /
+
+  DATA cg_aers_6bands/ &
+       ! bc soluble AS
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
+  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
+  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
+  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
+  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
+       ! pom soluble AS
+  0.687, 0.687, 0.687, 0.687, 0.687, 0.700, 0.710, 0.726, 0.736, 0.737, 0.738, 0.745, &
+  0.687, 0.687, 0.687, 0.687, 0.687, 0.700, 0.710, 0.726, 0.736, 0.737, 0.738, 0.745, &
+  0.658, 0.658, 0.658, 0.658, 0.658, 0.667, 0.674, 0.685, 0.692, 0.692, 0.693, 0.698, &
+  0.564, 0.564, 0.564, 0.564, 0.564, 0.566, 0.568, 0.571, 0.573, 0.573, 0.573, 0.574, &
+  0.363, 0.363, 0.363, 0.363, 0.363, 0.360, 0.357, 0.352, 0.350, 0.349, 0.349, 0.347, &
+  0.142, 0.142, 0.142, 0.142, 0.142, 0.139, 0.137, 0.133, 0.131, 0.131, 0.130, 0.129, &
+       ! sulfate AS 
+  0.675, 0.675, 0.675, 0.689, 0.701, 0.711, 0.720, 0.735, 0.748, 0.756, 0.760, 0.771, &
+  0.675, 0.675, 0.675, 0.689, 0.701, 0.711, 0.720, 0.735, 0.748, 0.756, 0.760, 0.771, &
+  0.653, 0.653, 0.653, 0.662, 0.670, 0.676, 0.683, 0.692, 0.701, 0.706, 0.709, 0.716, &
+  0.563, 0.563, 0.563, 0.565, 0.567, 0.569, 0.570, 0.573, 0.575, 0.576, 0.577, 0.579, &
+  0.362, 0.362, 0.362, 0.359, 0.356, 0.354, 0.352, 0.348, 0.345, 0.343, 0.342, 0.340, &
+  0.137, 0.137, 0.137, 0.135, 0.133, 0.132, 0.130, 0.128, 0.126, 0.125, 0.124, 0.122, &
+       ! sulfate coarse CS
+  0.803, 0.803, 0.803, 0.792, 0.783, 0.776, 0.769, 0.758, 0.747, 0.742, 0.738, 0.730, &
+  0.803, 0.803, 0.803, 0.792, 0.783, 0.776, 0.769, 0.758, 0.747, 0.742, 0.738, 0.730, &
+  0.799, 0.799, 0.799, 0.787, 0.777, 0.768, 0.760, 0.747, 0.736, 0.729, 0.725, 0.716, &
+  0.797, 0.797, 0.797, 0.782, 0.770, 0.760, 0.750, 0.735, 0.722, 0.714, 0.709, 0.698, &
+  0.810, 0.810, 0.810, 0.794, 0.781, 0.770, 0.759, 0.743, 0.728, 0.719, 0.714, 0.702, &
+  0.803, 0.803, 0.803, 0.790, 0.779, 0.770, 0.762, 0.748, 0.736, 0.729, 0.725, 0.715, &
+       ! seasalt seasalt Super Coarse Soluble (SS)
+  0.797, 0.799, 0.800, 0.803, 0.805, 0.823, 0.827, 0.830, 0.836, 0.838, 0.840, 0.845, &
+  0.788, 0.790, 0.792, 0.795, 0.798, 0.817, 0.821, 0.824, 0.830, 0.832, 0.835, 0.841, &
+  0.773, 0.776, 0.777, 0.782, 0.786, 0.805, 0.811, 0.814, 0.821, 0.824, 0.829, 0.835, &
+  0.746, 0.748, 0.752, 0.757, 0.764, 0.786, 0.792, 0.797, 0.806, 0.810, 0.815, 0.826, &
+  0.706, 0.710, 0.713, 0.718, 0.726, 0.752, 0.760, 0.766, 0.778, 0.783, 0.790, 0.806, &
+  0.681, 0.681, 0.681, 0.683, 0.685, 0.723, 0.728, 0.733, 0.740, 0.744, 0.751, 0.767, &
+       ! seasalt seasalt Coarse Soluble (CS)       
+  0.756, 0.759, 0.761, 0.766, 0.772, 0.792, 0.799, 0.802, 0.813, 0.816, 0.820, 0.828, &
+  0.736, 0.740, 0.742, 0.749, 0.756, 0.778, 0.785, 0.789, 0.800, 0.804, 0.810, 0.821, &
+  0.712, 0.714, 0.717, 0.724, 0.732, 0.757, 0.765, 0.770, 0.783, 0.788, 0.794, 0.809, &
+  0.690, 0.691, 0.693, 0.697, 0.703, 0.734, 0.742, 0.746, 0.758, 0.763, 0.769, 0.788, &
+  0.682, 0.682, 0.682, 0.683, 0.684, 0.725, 0.729, 0.733, 0.739, 0.742, 0.746, 0.759, &
+  0.669, 0.671, 0.672, 0.675, 0.678, 0.721, 0.727, 0.732, 0.737, 0.740, 0.743, 0.745, &
+       ! seasalt seasalt Accumulation Soluble (AS)
+  0.694, 0.694, 0.693, 0.690, 0.685, 0.733, 0.733, 0.735, 0.729, 0.728, 0.727, 0.732, &
+  0.685, 0.688, 0.690, 0.692, 0.693, 0.738, 0.742, 0.746, 0.745, 0.744, 0.742, 0.731, &
+  0.650, 0.656, 0.662, 0.672, 0.681, 0.722, 0.733, 0.740, 0.750, 0.753, 0.755, 0.748, &
+  0.561, 0.573, 0.584, 0.603, 0.626, 0.661, 0.681, 0.693, 0.721, 0.730, 0.741, 0.758, &
+  0.392, 0.408, 0.424, 0.453, 0.490, 0.517, 0.548, 0.567, 0.619, 0.637, 0.662, 0.716, &
+  0.144, 0.155, 0.167, 0.190, 0.223, 0.237, 0.268, 0.287, 0.353, 0.379, 0.418, 0.530  /
+
+  DATA cg_aeri_6bands/ &
+       ! dust insoluble CI
+  0.718, 0.718, 0.699, 0.661, 0.676, 0.670, &
+       ! bc insoluble AI
+  0.721, 0.721, 0.643, 0.513, 0.321, 0.153, &
+       ! pom insoluble AI
+  0.687, 0.687, 0.658, 0.564, 0.363, 0.142  /
+
+  DATA piz_aers_6bands/&
+       ! bc soluble AS
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
+  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
+  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
+  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
+  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
+       ! pom soluble AS
+  0.973, 0.973, 0.973, 0.973, 0.973, 0.977, 0.980, 0.984, 0.987, 0.988, 0.988, 0.990, &
+  0.973, 0.973, 0.973, 0.973, 0.973, 0.977, 0.980, 0.984, 0.987, 0.988, 0.988, 0.990, &
+  0.972, 0.972, 0.972, 0.972, 0.972, 0.976, 0.979, 0.984, 0.987, 0.987, 0.988, 0.990, &
+  0.940, 0.940, 0.940, 0.940, 0.940, 0.948, 0.955, 0.965, 0.972, 0.973, 0.973, 0.978, &
+  0.816, 0.816, 0.816, 0.816, 0.816, 0.839, 0.859, 0.888, 0.908, 0.910, 0.911, 0.925, &
+  0.663, 0.663, 0.663, 0.663, 0.663, 0.607, 0.562, 0.492, 0.446, 0.441, 0.437, 0.404, &
+       ! sulfate 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, &
+  1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  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.988, 0.988, 0.988, 0.989, 0.990, 0.990, 0.991, 0.992, 0.993, 0.993, 0.994, 0.994, &
+  0.256, 0.256, 0.256, 0.263, 0.268, 0.273, 0.277, 0.284, 0.290, 0.294, 0.296, 0.301, &
+       ! sulfate coarse 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, &
+  1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  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.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, &
+  0.877, 0.877, 0.877, 0.873, 0.870, 0.867, 0.864, 0.860, 0.856, 0.854, 0.852, 0.849, &
+       ! seasalt seasalt Super Coarse Soluble (SS)
+  0.983, 0.982, 0.982, 0.982, 0.981, 0.992, 0.993, 0.994, 0.995, 0.996, 0.997, 0.998, &
+  0.984, 0.983, 0.983, 0.983, 0.982, 0.992, 0.993, 0.995, 0.996, 0.996, 0.997, 0.998, &
+  0.986, 0.985, 0.985, 0.985, 0.984, 0.993, 0.994, 0.995, 0.996, 0.997, 0.997, 0.998, &
+  0.989, 0.989, 0.988, 0.988, 0.987, 0.995, 0.995, 0.996, 0.997, 0.997, 0.998, 0.999, &
+  0.994, 0.993, 0.993, 0.992, 0.992, 0.997, 0.997, 0.997, 0.998, 0.998, 0.998, 0.999, &
+  0.997, 0.997, 0.997, 0.997, 0.996, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, &
+       ! seasalt seasalt Coarse Soluble (CS)      
+  0.988, 0.988, 0.988, 0.987, 0.986, 0.994, 0.995, 0.996, 0.996, 0.997, 0.998, 0.998, &
+  0.990, 0.990, 0.990, 0.989, 0.988, 0.995, 0.996, 0.996, 0.997, 0.997, 0.998, 0.999, &
+  0.993, 0.993, 0.992, 0.992, 0.991, 0.996, 0.997, 0.997, 0.998, 0.998, 0.998, 0.999, &
+  0.996, 0.995, 0.995, 0.994, 0.994, 0.998, 0.998, 0.998, 0.998, 0.999, 0.999, 0.999, &
+  0.998, 0.997, 0.997, 0.997, 0.997, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, &
+  0.999, 0.999, 0.998, 0.998, 0.998, 0.999, 0.999, 1.000, 1.000, 1.000, 1.000, 1.000, &
+       ! seasalt seasalt Accumulation Soluble (AS)
+  0.999, 0.999, 0.998, 0.998, 0.998, 0.999, 0.999, 0.999, 0.999, 1.000, 1.000, 1.000, &
+  0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  0.999, 0.999, 0.999, 0.999, 0.999, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  0.999, 0.999, 0.999, 0.999, 0.999, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  0.997, 0.997, 0.998, 0.998, 0.998, 0.999, 0.999, 1.000, 1.000, 1.000, 1.000, 1.000, &
+  0.985, 0.987, 0.988, 0.991, 0.993, 0.997, 0.998, 0.998, 0.999, 0.999, 1.000, 1.000  /
+
+  DATA piz_aeri_6bands/ &
+       ! dust insoluble CI
+  0.954, 0.954, 0.965, 0.981, 0.996, 0.990, &
+       ! bc insoluble AI
+  0.460, 0.460, 0.445, 0.394, 0.267, 0.121, &
+       ! pom insoluble AI
+  0.973, 0.973, 0.972, 0.940, 0.816, 0.663 /
+
+!----BEGINNING OF CALCULATIONS
+
+  spsol = 0
+  spinsol = 0 
+  if(NSW.NE.nbands_rrtm) then
+     print *,'Erreur NSW doit etre egal a 6 pour cette routine'
+     stop
+  endif
+
+  DO k=1, klev
+    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
+!CDIR UNROLL=naero_spc
+      mass_temp_pi(i,k,:) = m_allaer_pi(i,k,:) / zrho / 1.e+9
+      zdp1(i,k)=pdel(i,k)/(RG*delt)      ! air mass auxiliary  variable --> zdp1 [kg/(m^2 *s)]
+    ENDDO
+  ENDDO
+
+  IF (flag_aerosol .EQ. 1) THEN 
+     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
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASBCM
+     aerosol_name(2) = id_AIBCM
+  ELSEIF (flag_aerosol .EQ. 3) THEN 
+     nb_aer = 2
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASPOMM
+     aerosol_name(2) = id_AIPOMM
+  ELSEIF (flag_aerosol .EQ. 4) THEN 
+     nb_aer = 3
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_CSSSM
+     aerosol_name(2) = id_SSSSM
+     aerosol_name(3) = id_ASSSM
+  ELSEIF (flag_aerosol .EQ. 5) THEN 
+     nb_aer = 1
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_CIDUSTM
+  ELSEIF (flag_aerosol .EQ. 6) THEN 
+     nb_aer = 10
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASSO4M      
+     aerosol_name(2) = id_ASBCM
+     aerosol_name(3) = id_AIBCM
+     aerosol_name(4) = id_ASPOMM
+     aerosol_name(5) = id_AIPOMM
+     aerosol_name(6) = id_CSSSM
+     aerosol_name(7) = id_SSSSM
+     aerosol_name(8) = id_ASSSM
+     aerosol_name(9) = id_CIDUSTM
+     aerosol_name(10)= id_CSSO4M
+  ENDIF
+
+  !
+  ! loop over modes, use of precalculated nmd and corresponding sigma
+  !    loop over wavelengths
+  !    for each mass species in mode
+  !      interpolate from Sext to retrieve Sext_at_gridpoint_per_species
+  !      compute optical_thickness_at_gridpoint_per_species
+
+!!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(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
+     IF (aerosol_name(m).EQ.id_ASBCM) THEN
+        soluble=.TRUE.
+        spsol=1
+     ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN 
+        soluble=.TRUE.
+        spsol=2 
+     ELSEIF (aerosol_name(m).EQ.id_ASSO4M) THEN 
+        soluble=.TRUE.
+        spsol=3
+        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
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+     ELSEIF (aerosol_name(m).EQ.id_SSSSM) THEN 
+         soluble=.TRUE.
+         spsol=5
+     ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN 
+         soluble=.TRUE.
+         spsol=6
+     ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN
+         soluble=.TRUE.
+         spsol=7
+     ELSEIF (aerosol_name(m).EQ.id_CIDUSTM) THEN 
+         soluble=.FALSE.
+         spinsol=1
+     ELSEIF  (aerosol_name(m).EQ.id_AIBCM) THEN 
+         soluble=.FALSE.
+         spinsol=2
+     ELSEIF (aerosol_name(m).EQ.id_AIPOMM) THEN 
+         soluble=.FALSE.
+         spinsol=3
+     ELSE 
+         CYCLE
+     ENDIF
+
+    id=aerosol_name(m)
+    used_aer(id)=.TRUE.
+
+    IF (soluble) THEN
+
+       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
+
+           DO inu=1,NSW
+
+             tau_ae2b_int= alpha_aers_6bands(RH_num(i,k),inu,spsol)+ & 
+                           delta(i,k)* (alpha_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
+                           alpha_aers_6bands(RH_num(i,k),inu,spsol))
+                   
+             piz_ae2b_int = piz_aers_6bands(RH_num(i,k),inu,spsol) + & 
+                            delta(i,k)* (piz_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
+                            piz_aers_6bands(RH_num(i,k),inu,spsol))
+                   
+             cg_ae2b_int = cg_aers_6bands(RH_num(i,k),inu,spsol) + & 
+                           delta(i,k)* (cg_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 
+                           cg_aers_6bands(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
+        
+     ELSE    ! For all aerosol insoluble components
+
+       DO k=1, KLEV
+         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
+
+           DO inu=1,NSW
+             tau_ae2b_int = alpha_aeri_6bands(inu,spinsol)
+             piz_ae2b_int = piz_aeri_6bands(inu,spinsol)
+             cg_ae2b_int = cg_aeri_6bands(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 / insoluble
+
+  ENDDO  ! nb_aer  
+
+  DO m=1,naero_tot
+    IF (.NOT. used_aer(m)) THEN
+      tau_ae(:,:,m,:)=0.
+      tau_ae_pi(:,:,m,:)=0.
+      piz_ae(:,:,m,:)=0.
+      cg_ae(:,:,m,:)=0.
+    ENDIF
+  ENDDO
+
+  DO inu=1, NSW 
+     DO k=1, KLEV
+       DO i=1, KLON
+!--anthropogenic aerosol
+         tau_allaer(i,k,2,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,2,inu)=MAX(tau_allaer(i,k,2,inu),1e-5)
+
+         piz_allaer(i,k,2,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,2,inu)
+         piz_allaer(i,k,2,inu)=MAX(piz_allaer(i,k,2,inu),0.1)
+
+         cg_allaer(i,k,2,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)+ &
+                               tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,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_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_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu)*cg_ae(i,k,id_CIDUSTM,inu))/ &
+                               (tau_allaer(i,k,2,inu)*piz_allaer(i,k,2,inu))
+
+!--natural aerosol
+         tau_allaer(i,k,1,inu)=tau_ae_pi(i,k,id_ASSO4M,inu)+tau_ae_pi(i,k,id_CSSO4M,inu)+ &
+                               tau_ae_pi(i,k,id_ASBCM,inu)+tau_ae_pi(i,k,id_AIBCM,inu)+   &
+                               tau_ae_pi(i,k,id_ASPOMM,inu)+tau_ae_pi(i,k,id_AIPOMM,inu)+ &
+                               tau_ae_pi(i,k,id_ASSSM,inu)+tau_ae_pi(i,k,id_CSSSM,inu)+   &
+                               tau_ae_pi(i,k,id_SSSSM,inu)+ tau_ae_pi(i,k,id_CIDUSTM,inu)
+         tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),1e-5)
+
+         piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)+ &
+                                tau_ae_pi(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)+ &
+                                tau_ae_pi(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)+ &
+                                tau_ae_pi(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)+ &
+                                tau_ae_pi(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)+ &
+                                tau_ae_pi(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)+ &
+                                tau_ae_pi(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)+ &
+                                tau_ae_pi(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)+ &
+                                tau_ae_pi(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)+ &
+                                tau_ae_pi(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu)) &
+                                /tau_allaer(i,k,1,inu)
+         piz_allaer(i,k,1,inu)=MAX(piz_allaer(i,k,1,inu),0.1)
+
+         cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ &
+                               tau_ae_pi(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)*cg_ae(i,k,id_CSSO4M,inu)+ &
+                               tau_ae_pi(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)*cg_ae(i,k,id_ASBCM,inu)+ &
+                               tau_ae_pi(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu)+ &
+                               tau_ae_pi(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)*cg_ae(i,k,id_ASPOMM,inu)+ &
+                               tau_ae_pi(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu)+ &
+                               tau_ae_pi(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)*cg_ae(i,k,id_ASSSM,inu)+ &
+                               tau_ae_pi(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu)+ &
+                               tau_ae_pi(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu)+ &
+                               tau_ae_pi(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu)*cg_ae(i,k,id_CIDUSTM,inu))/ &
+                               (tau_allaer(i,k,1,inu)*piz_allaer(i,k,1,inu))
+
+        ENDDO
+      ENDDO
+    ENDDO
+   
+!--???????
+  inu=1
+  DO i=1, KLON
+     absvisaer(i)=SUM((1-piz_allaer(i,:,:,inu))*tau_allaer(i,:,:,inu))
+  END DO	
+
+  DEALLOCATE(aerosol_name) 
+
+END SUBROUTINE AEROPT_6BANDS_RRTM
Index: LMDZ5/trunk/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90	(revision 2004)
+++ LMDZ5/trunk/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90	(revision 2004)
@@ -0,0 +1,215 @@
+! $Id$
+!
+SUBROUTINE readaerosol_optic_rrtm(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, &
+     tausum_aero, tau3d_aero )
+
+! This routine will :
+! 1) recevie the aerosols(already read and interpolated) corresponding to flag_aerosol
+! 2) calculate the optical properties for the aerosols
+!
+  
+  USE dimphy
+  USE aero_mod
+  USE phys_local_var_mod, only: sconcso4,sconcoa,sconcbc,sconcss,sconcdust, &
+      concso4,concoa,concbc,concss,concdust,loadso4,loadoa,loadbc,loadss,loaddust, &
+      load_tmp1,load_tmp2,load_tmp3,load_tmp4,load_tmp5,load_tmp6,load_tmp7
+  IMPLICIT NONE
+  include "clesphys.h"
+
+! Input arguments
+!****************************************************************************************
+  LOGICAL, INTENT(IN)                      :: debut
+  LOGICAL, INTENT(IN)                      :: new_aod
+  INTEGER, INTENT(IN)                      :: flag_aerosol
+  INTEGER, INTENT(IN)                      :: itap
+  REAL, INTENT(IN)                         :: rjourvrai
+  REAL, INTENT(IN)                         :: pdtphys
+  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
+  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
+  REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_seri
+  REAL, DIMENSION(klon,klev), INTENT(IN)   :: rhcl   ! humidite relative ciel clair
+  REAL, DIMENSION(klev), INTENT(IN)        :: presnivs
+
+! Output arguments
+!****************************************************************************************
+  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero    ! Total mass for all soluble aerosols
+  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero_pi !     -"-     preindustrial values
+  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: tau_aero    ! Aerosol optical thickness
+  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: piz_aero    ! Single scattering albedo aerosol
+  REAL, DIMENSION(klon,klev,2,NSW), INTENT(OUT) :: cg_aero     ! asymmetry parameter aerosol
+!  REAL, DIMENSION(klon,nwave,naero_spc), INTENT(OUT)       :: tausum_aero
+!  REAL, DIMENSION(klon,klev,nwave,naero_spc), INTENT(OUT)  :: tau3d_aero
+!--correction minibug OB
+  REAL, DIMENSION(klon,nwave,naero_tot), INTENT(OUT)       :: tausum_aero
+  REAL, DIMENSION(klon,klev,nwave,naero_tot), INTENT(OUT)  :: tau3d_aero
+
+! Local variables
+!****************************************************************************************
+  REAL, DIMENSION(klon)        :: aerindex ! POLDER aerosol index 
+  REAL, DIMENSION(klon,klev)   :: sulfate  ! SO4 aerosol concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: bcsol    ! BC soluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: bcins    ! BC insoluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: pomsol   ! POM soluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: pomins   ! POM insoluble concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: cidust    ! DUST aerosol concentration  [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: sscoarse  ! SS Coarse concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: sssupco   ! SS Super Coarse concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: ssacu     ! SS Acumulation concentration [ug/m3]
+  REAL, DIMENSION(klon,klev)   :: sulfate_pi
+  REAL, DIMENSION(klon,klev)   :: bcsol_pi
+  REAL, DIMENSION(klon,klev)   :: bcins_pi
+  REAL, DIMENSION(klon,klev)   :: pomsol_pi
+  REAL, DIMENSION(klon,klev)   :: pomins_pi
+  REAL, DIMENSION(klon,klev)   :: cidust_pi
+  REAL, DIMENSION(klon,klev)   :: sscoarse_pi
+  REAL, DIMENSION(klon,klev)   :: sssupco_pi
+  REAL, DIMENSION(klon,klev)   :: ssacu_pi
+  REAL, DIMENSION(klon,klev)   :: pdel
+  REAL, DIMENSION(klon,klev,naero_spc) :: m_allaer
+  REAL, DIMENSION(klon,klev,naero_spc) :: m_allaer_pi !RAF  
+!  REAL, DIMENSION(klon,naero_tot)      :: fractnat_allaer !RAF delete??
+
+  INTEGER :: k, i
+  
+!****************************************************************************************
+! 1) Get aerosol mass
+!    
+!****************************************************************************************
+! Read and interpolate sulfate
+  IF ( flag_aerosol .EQ. 1 .OR. &
+       flag_aerosol .EQ. 6 ) THEN 
+
+     CALL readaerosol_interp(id_ASSO4M, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sulfate, sulfate_pi,loadso4)
+  ELSE
+     sulfate(:,:) = 0. ; sulfate_pi(:,:) = 0.
+     loadso4=0.
+  END IF
+
+! Read and interpolate bcsol and bcins
+  IF ( flag_aerosol .EQ. 2 .OR. &
+       flag_aerosol .EQ. 6 ) THEN 
+
+     ! Get bc aerosol distribution 
+     CALL readaerosol_interp(id_ASBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi, load_tmp1 )
+     CALL readaerosol_interp(id_AIBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi, load_tmp2 )
+     loadbc(:)=load_tmp1(:)+load_tmp2(:)
+  ELSE
+     bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0.
+     bcins(:,:) = 0. ; bcins_pi(:,:) = 0.
+     loadbc=0.
+  END IF
+
+
+! Read and interpolate pomsol and pomins
+  IF ( flag_aerosol .EQ. 3 .OR. &
+       flag_aerosol .EQ. 6 ) THEN
+
+     CALL readaerosol_interp(id_ASPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp3)
+     CALL readaerosol_interp(id_AIPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp4)
+     loadoa(:)=load_tmp3(:)+load_tmp4(:)
+  ELSE
+     pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
+     pomins(:,:) = 0. ; pomins_pi(:,:) = 0.
+     loadoa=0.
+  END IF
+
+
+! Read and interpolate csssm, ssssm, assssm
+  IF (flag_aerosol .EQ. 4 .OR. &
+      flag_aerosol .EQ. 6 ) THEN 
+
+      CALL readaerosol_interp(id_SSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi, load_tmp5) 
+      CALL readaerosol_interp(id_CSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp6) 
+      CALL readaerosol_interp(id_ASSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, ssacu, ssacu_pi, load_tmp7) 
+     loadss(:)=load_tmp5(:)+load_tmp6(:)+load_tmp7(:)
+  ELSE
+     sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0. 
+     ssacu(:,:)    = 0. ; ssacu_pi(:,:) = 0. 
+     sssupco(:,:)  = 0. ; sssupco_pi = 0. 
+     loadss=0.
+  ENDIF
+
+! Read and interpolate cidustm
+  IF (flag_aerosol .EQ. 5 .OR.  &
+      flag_aerosol .EQ. 6 ) THEN 
+
+      CALL readaerosol_interp(id_CIDUSTM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi, loaddust) 
+
+  ELSE
+      cidust(:,:) = 0. ; cidust_pi(:,:) = 0. 
+      loaddust=0.
+  ENDIF
+
+!
+! Store all aerosols in one variable
+!
+  m_allaer(:,:,id_ASBCM)  = bcsol(:,:)        ! ASBCM
+  m_allaer(:,:,id_ASPOMM) = pomsol(:,:)       ! ASPOMM
+  m_allaer(:,:,id_ASSO4M) = sulfate(:,:)      ! ASSO4M (= SO4) 
+  m_allaer(:,:,id_CSSO4M) = 0.                ! CSSO4M 
+  m_allaer(:,:,id_SSSSM)  = sssupco(:,:)      ! SSSSM
+  m_allaer(:,:,id_CSSSM)  = sscoarse(:,:)     ! CSSSM
+  m_allaer(:,:,id_ASSSM)  = ssacu(:,:)        ! ASSSM
+  m_allaer(:,:,id_CIDUSTM)= cidust(:,:)       ! CIDUSTM
+  m_allaer(:,:,id_AIBCM)  = bcins(:,:)        ! AIBCM
+  m_allaer(:,:,id_AIPOMM) = pomins(:,:)       ! AIPOMM
+
+!RAF
+  m_allaer_pi(:,:,1)  = bcsol_pi(:,:)        ! ASBCM pre-ind
+  m_allaer_pi(:,:,2)  = pomsol_pi(:,:)       ! ASPOMM pre-ind
+  m_allaer_pi(:,:,3)  = sulfate_pi(:,:)      ! ASSO4M (= SO4) pre-ind
+  m_allaer_pi(:,:,4)  = 0.                   ! CSSO4M pre-ind
+  m_allaer_pi(:,:,5)  = sssupco_pi(:,:)      ! SSSSM pre-ind
+  m_allaer_pi(:,:,6)  = sscoarse_pi(:,:)     ! CSSSM pre-ind
+  m_allaer_pi(:,:,7)  = ssacu_pi(:,:)        ! ASSSM pre-ind
+  m_allaer_pi(:,:,8)  = cidust_pi(:,:)       ! CIDUSTM pre-ind
+  m_allaer_pi(:,:,9)  = bcins_pi(:,:)        ! AIBCM pre-ind
+  m_allaer_pi(:,:,10) = pomins_pi(:,:)       ! AIPOMM pre-ind
+
+!
+! Calculate the total mass of all soluble aersosols
+!
+  mass_solu_aero(:,:)    = sulfate(:,:)    + bcsol(:,:)    + pomsol(:,:) !   + &
+!       sscoarse(:,:)    + ssacu(:,:)    + sssupco(:,:) 
+  mass_solu_aero_pi(:,:) = sulfate_pi(:,:) + bcsol_pi(:,:) + pomsol_pi(:,:) ! + &
+!       sscoarse_pi(:,:) + ssacu_pi(:,:) + sssupco_pi(:,:)
+
+!****************************************************************************************
+! 2) Calculate optical properties for the aerosols
+!
+!****************************************************************************************
+  DO k = 1, klev
+     DO i = 1, klon
+        pdel(i,k) = paprs(i,k) - paprs (i,k+1)
+     END DO
+  END DO
+
+  CALL aeropt_6bands_rrtm( &
+          pdel, m_allaer, pdtphys, rhcl, & 
+          tau_aero, piz_aero, cg_aero,   &
+          m_allaer_pi, flag_aerosol, &
+          pplay, t_seri, presnivs) 
+     
+     ! aeropt_5wv only for validation and diagnostics.
+     CALL aeropt_5wv_rrtm(                    &
+          pdel, m_allaer,                &
+          pdtphys, rhcl, aerindex,       & 
+          flag_aerosol, pplay, t_seri,   &
+          tausum_aero, tau3d_aero )
+
+! Diagnostics calculation for CMIP5 protocol
+  sconcso4(:)=m_allaer(:,1,id_ASSO4M)*1.e-9
+  sconcoa(:)=(m_allaer(:,1,id_ASPOMM)+m_allaer(:,1,id_AIPOMM))*1.e-9
+  sconcbc(:)=(m_allaer(:,1,id_ASBCM)+m_allaer(:,1,id_AIBCM))*1.e-9
+  sconcss(:)=(m_allaer(:,1,id_ASSSM)+m_allaer(:,1,id_CSSSM)+m_allaer(:,1,id_SSSSM))*1.e-9
+  sconcdust(:)=m_allaer(:,1,id_CIDUSTM)*1.e-9
+  concso4(:,:)=m_allaer(:,:,id_ASSO4M)*1.e-9
+  concoa(:,:)=(m_allaer(:,:,id_ASPOMM)+m_allaer(:,:,id_AIPOMM))*1.e-9
+  concbc(:,:)=(m_allaer(:,:,id_ASBCM)+m_allaer(:,:,id_AIBCM))*1.e-9
+  concss(:,:)=(m_allaer(:,:,id_ASSSM)+m_allaer(:,:,id_CSSSM)+m_allaer(:,:,id_SSSSM))*1.e-9
+  concdust(:,:)=m_allaer(:,:,id_CIDUSTM)*1.e-9
+
+END SUBROUTINE readaerosol_optic_rrtm
Index: LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato_rrtm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato_rrtm.F90	(revision 2004)
+++ LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato_rrtm.F90	(revision 2004)
@@ -0,0 +1,172 @@
+subroutine readaerosolstrato_rrtm(debut)
+
+    use netcdf95, only: nf95_close, nf95_gw_var, nf95_inq_dimid, & 
+                        nf95_inq_varid, nf95_open
+    use netcdf, only: nf90_get_var, nf90_noerr, nf90_nowrite
+
+    USE phys_cal_mod, ONLY : mth_cur
+    USE mod_grid_phy_lmdz
+    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+    USE mod_phys_lmdz_para 
+    USE phys_state_var_mod
+    USE phys_local_var_mod
+    USE aero_mod
+    USE dimphy
+
+    implicit none
+
+    include "YOMCST.h"
+    include "dimensions.h"
+
+! Variable input
+    logical debut
+
+! Variables locales
+    integer n_lat   ! number of latitudes in the input data
+    integer n_lon   ! number of longitudes in the input data
+    integer n_lev   ! number of levels in the input data
+    integer n_month ! number of months in the input data
+    real, pointer:: latitude(:)
+    real, pointer:: longitude(:)
+    real, pointer:: time(:)
+    real, pointer:: lev(:)
+    integer k, band, wave
+    integer, save :: mth_pre
+
+    real, allocatable, dimension(:,:), save :: tau_aer_strat
+!$OMP THREADPRIVATE(tau_aer_strat)
+
+! Champs reconstitues
+    real, allocatable:: tauaerstrat(:, :, :, :)
+    real, allocatable:: tauaerstrat_mois(:, :, :)
+    real, allocatable:: tauaerstrat_mois_glo(:, :)
+    real, allocatable:: tauaerstrat_mois_glo_bands(:,:,:)
+
+! For NetCDF:
+    integer ncid_in  ! IDs for input files
+    integer varid, ncerr
+
+! Stratospheric aerosols optical properties
+! alpha_strat over the 2 bands is normalised by the 550 nm extinction coefficient
+! alpha_strat_wave is *not* normalised by the 550 nm extinction coefficient
+    real, dimension(nbands_rrtm) :: alpha_strat, piz_strat, cg_strat
+    data alpha_strat/0.938538969, 0.990073204, 0.992904723, 0.829215884, 0.439313501, 0.156857833/
+    data cg_strat   /0.699142992, 0.716326416, 0.735462785, 0.736726701, 0.712068975, 0.575097859/
+    data piz_strat  /1.000000000, 1.000000000, 1.000000000, 1.000000000, 0.997781098, 0.452584684/
+    real, dimension(nwave) :: alpha_strat_wave
+    data alpha_strat_wave/3.36780953,3.34667683,3.20444202,3.0293026,2.82108808/
+
+!--------------------------------------------------------
+
+    IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev))
+
+    IF (is_mpi_root) THEN
+
+    IF (debut.OR.mth_cur.NE.mth_pre) THEN
+
+    IF (nbands_rrtm.NE.6) THEN 
+        print *,'nbands_rrtm doit etre egal a 6 dans readaerosolstrat_rrtm'
+        STOP
+    ENDIF
+
+    CALL nf95_open("taustrat.nc", nf90_nowrite, ncid_in)
+
+    CALL nf95_inq_varid(ncid_in, "LEV", varid)
+    CALL nf95_gw_var(ncid_in, varid, lev)
+    n_lev = size(lev)
+    IF (n_lev.NE.klev) THEN 
+       print *,'Le nombre de niveaux n est pas egal a klev'
+       STOP
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "LAT", varid)
+    CALL nf95_gw_var(ncid_in, varid, latitude)
+    n_lat = size(latitude)
+    print *, 'LAT aerosol strato=', n_lat, latitude
+    IF (n_lat.NE.jjm+1) THEN 
+       print *,'Le nombre de lat n est pas egal a jjm+1'
+       STOP
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "LON", varid)
+    CALL nf95_gw_var(ncid_in, varid, longitude)
+    n_lon = size(longitude)
+    print *, 'LON aerosol strato=', n_lon, longitude
+    IF (n_lon.NE.iim) THEN 
+       print *,'Le nombre de lon n est pas egal a iim'
+       STOP
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "TIME", varid)
+    CALL nf95_gw_var(ncid_in, varid, time)
+    n_month = size(time)
+    print *, 'TIME aerosol strato=', n_month, time
+    IF (n_month.NE.12) THEN 
+       print *,'Le nombre de month n est pas egal a 12'
+       STOP
+    ENDIF
+
+    ALLOCATE(tauaerstrat(n_lon, n_lat, n_lev, n_month))
+    ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev))
+    ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev))
+    ALLOCATE(tauaerstrat_mois_glo_bands(klon_glo, n_lev,nbands_rrtm))
+
+!--reading stratospheric AOD at 550 nm
+    CALL nf95_inq_varid(ncid_in, "TAUSTRAT", varid)
+    ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
+    print *,'code erreur readaerosolstrato=', ncerr, varid
+
+    CALL nf95_close(ncid_in)
+
+!---select the correct month
+    IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
+      print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
+    ENDIF
+    tauaerstrat_mois(:,:,:) = tauaerstrat(:,:,:,mth_cur)
+
+!---reduce to a klon_glo grid 
+    CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
+
+!--scatter on all proc
+    CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
+
+    DEALLOCATE(tauaerstrat)
+    DEALLOCATE(tauaerstrat_mois)
+    DEALLOCATE(tauaerstrat_mois_glo)
+  
+    mth_pre=mth_cur
+
+    ENDIF !--debut ou nouveau mois
+
+    ENDIF !--is_mpi_root
+
+!--total vertical aod at the 5 wavelengths
+    DO wave=1, nwave
+    DO k=1, klev
+    tausum_aero(:,wave,id_strat)=tausum_aero(:,wave,id_strat)+tau_aer_strat(:,k)*alpha_strat_wave(wave)/alpha_strat_wave(2)
+    ENDDO
+    ENDDO
+
+!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
+    DO band=1, nbands_rrtm
+!--anthropogenic aerosols bands 1 to nbands_rrtm 
+    cg_aero_rrtm(:,:,2,band)  = ( cg_aero_rrtm(:,:,2,band)*piz_aero_rrtm(:,:,2,band)*tau_aero_rrtm(:,:,2,band) + &
+                                  cg_strat(band)*piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /        &
+                             MAX( piz_aero_rrtm(:,:,2,band)*tau_aero_rrtm(:,:,2,band) +                          &
+                                  piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:), 1.e-15 )
+    piz_aero_rrtm(:,:,2,band)  = ( piz_aero_rrtm(:,:,2,band)*tau_aero_rrtm(:,:,2,band) +                         &
+                              piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /                           &
+                              MAX( tau_aero_rrtm(:,:,2,band) + alpha_strat(band)*tau_aer_strat(:,:), 1.e-15 )
+    tau_aero_rrtm(:,:,2,band)  = tau_aero_rrtm(:,:,2,band) + alpha_strat(band)*tau_aer_strat(:,:)
+!--natural aerosols bands 1 to nbands_rrtm
+    cg_aero_rrtm(:,:,1,band)  = ( cg_aero_rrtm(:,:,1,band)*piz_aero_rrtm(:,:,1,band)*tau_aero_rrtm(:,:,1,band) + &
+                             cg_strat(band)*piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /             &
+                             MAX( piz_aero_rrtm(:,:,1,band)*tau_aero_rrtm(:,:,1,band) +                          &
+                                  piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:), 1.e-15 )
+    piz_aero_rrtm(:,:,1,band)  = ( piz_aero_rrtm(:,:,1,band)*tau_aero_rrtm(:,:,1,band) +                         & 
+                              piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /                           &
+                              MAX( tau_aero_rrtm(:,:,1,band) + alpha_strat(band)*tau_aer_strat(:,:),1.e-15 )
+    tau_aero_rrtm(:,:,1,band)  = tau_aero_rrtm(:,:,1,band) + alpha_strat(band)*tau_aer_strat(:,:)
+    ENDDO
+
+end subroutine readaerosolstrato_rrtm
Index: LMDZ5/trunk/libf/phylmd/rrtm/recmwf_aero.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/rrtm/recmwf_aero.F90	(revision 2004)
+++ LMDZ5/trunk/libf/phylmd/rrtm/recmwf_aero.F90	(revision 2004)
@@ -0,0 +1,792 @@
+!OPTIONS XOPT(NOEVAL)
+SUBROUTINE RECMWF_AERO (KST, KEND, KPROMA, KTDIA , KLEV,&
+ & KMODE,&
+ & PALBD , PALBP , PAPRS , PAPRSF , PCCO2 , PCLFR,&
+ & PQO3  , PAER  , PDP   , PEMIS  , PMU0,&
+ & PQ    , PQS   , PQIWP , PQLWP , PSLM   , PT    , PTS,&
+ & PREF_LIQ, PREF_ICE,&
+!--OB
+ & PREF_LIQ_PI, PREF_ICE_PI,&
+!--fin
+ & PEMTD , PEMTU , PTRSO,&
+ & PTH   , PCTRSO, PCEMTR, PTRSOD,&
+ & PLWFC, PLWFT, PSWFC, PSWFT, PSFSWDIR, PSFSWDIF,&
+ & PFSDNN, PFSDNV,&  
+ & PPIZA_TOT,PCGA_TOT,PTAU_TOT, & 
+!--OB
+ & PPIZA_NAT,PCGA_NAT,PTAU_NAT, & 
+!--fin OB 
+ & PFLUX,PFLUC,&
+ & PFSDN ,PFSUP , PFSCDN , PFSCUP,&
+!--OB diagnostics
+ & PTOPSWADAERO,PSOLSWADAERO,&
+ & PTOPSWAD0AERO,PSOLSWAD0AERO,&
+ & PTOPSWAIAERO,PSOLSWAIAERO,&
+ & PTOPSWCFAERO,PSOLSWCFAERO,&
+ & ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
+!--fin
+
+!**** *RECMWF* - METEO-FRANCE RADIATION INTERFACE TO ECMWF RADIATION SCHEME
+
+!     PURPOSE.
+!     --------
+!           SIMPLE INTERFACE TO RADLSW (NO INTERPOLATION)
+
+!**   INTERFACE.
+!     ----------
+
+!     EXPLICIT ARGUMENTS :
+!        --------------------
+! KST    : START INDEX OF DATA IN KPROMA-LONG VECTOR
+! KEND   : END   INDEX OF DATA IN KPROMA-LONG VECTOR
+! KPROMA : VECTOR LENGTH
+! KTDIA  : INDEX OF TOP LEVEL FROM WHICH COMPUTATIONS ARE ACTIVE
+! KLEV   : NUMBER OF LEVELS
+! PAER   : (KPROMA,KLEV ,6)     ; OPTICAL THICKNESS OF THE AEROSOLS
+! PALBD  : (KPROMA,NSW)         ; DIFFUSE ALBEDO IN THE 2 SW INTERVALS
+! PALBP  : (KPROMA,NSW)         ; PARALLEL ALBEDO IN THE 2 SW INTERVALS
+! PAPRS  : (KPROMA,KLEV+1)      ; HALF LEVEL PRESSURE
+! PAPRSF : (KPROMA,KLEV )       ; FULL LEVEL PRESSURE
+! PCCO2  :                      ; CONCENTRATION IN CO2 (PA/PA)
+! PCLFR  : (KPROMA,KLEV )       ; CLOUD FRACTIONAL COVER
+! PQO3   : (KPROMA,KLEV )       ; OZONE MIXING RATIO (MASS)
+! PDP    : (KPROMA,KLEV)        ; LAYER PRESSURE THICKNESS
+! PEMIS  : (KPROMA)             ; SURFACE EMISSIVITY
+! PMU0   : (KPROMA)             ; SOLAR ANGLE
+! PQ     : (KPROMA,KLEV )       ; SPECIFIC HUMIDITY PA/PA
+! PQS    : (KPROMA,KLEV )       ; SATURATION SPECIFIC HUMIDITY PA/PA
+! PQIWP  : (KPROMA,KLEV )       ; ICE    WATER KG/KG
+! PQLWP  : (KPROMA,KLEV )       ; LIQUID WATER KG/KG
+! PSLM   : (KPROMA)             ; LAND-SEA MASK
+! PT     : (KPROMA,KLEV)        ; FULL LEVEL TEMPERATURE
+! PTS    : (KPROMA)             ; SURFACE TEMPERATURE
+! PPIZA_TOT  : (KPROMA,KLEV,NSW); Single scattering albedo of total aerosol
+! PCGA_TOT   : (KPROMA,KLEV,NSW); Assymetry factor for total aerosol
+! PTAU_TOT: (KPROMA,KLEV,NSW)   ; Optical depth of total aerosol
+! PREF_LIQ (KPROMA,KLEV)        ; Liquid droplet radius (um) - present-day
+! PREF_ICE (KPROMA,KLEV)        ; Ice crystal radius (um) - present-day
+!--OB
+! PREF_LIQ_PI (KPROMA,KLEV)     ; Liquid droplet radius (um) - pre-industrial 
+! PREF_ICE_PI (KPROMA,KLEV)     ; Ice crystal radius (um) - pre-industrial 
+! ok_ade---input-L- apply the Aerosol Direct Effect or not?
+! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
+! flag_aerosol-input-I- aerosol flag from 0 to 6
+! flag_aerosol_strat-input-I- use stratospheric aerosols flag (T/F)
+! PPIZA_NAT  : (KPROMA,KLEV,NSW); Single scattering albedo of natural aerosol
+! PCGA_NAT   : (KPROMA,KLEV,NSW); Assymetry factor for natural aerosol
+! PTAU_NAT: (KPROMA,KLEV,NSW)   ; Optical depth of natural aerosol
+!--fin OB
+
+!     ==== OUTPUTS ===
+! PEMTD (KPROMA,KLEV+1)         ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
+! PEMTU (KPROMA,KLEV+1)         ; TOTAL UPWARD   LONGWAVE EMISSIVITY
+! PTRSO (KPROMA,KLEV+1)         ; TOTAL SHORTWAVE TRANSMISSIVITY
+! PTH   (KPROMA,KLEV+1)         ; HALF LEVEL TEMPERATURE
+! PCTRSO(KPROMA,2)              ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
+! PCEMTR(KPROMA,2)              ; CLEAR-SKY NET LONGWAVE EMISSIVITY
+! PTRSOD(KPROMA)                ; TOTAL-SKY SURFACE SW TRANSMISSITY
+! PLWFC (KPROMA,2)              ; CLEAR-SKY LONGWAVE FLUXES
+! PLWFT (KPROMA,KLEV+1)         ; TOTAL-SKY LONGWAVE FLUXES
+! PSWFC (KPROMA,2)              ; CLEAR-SKY SHORTWAVE FLUXES
+! PSWFT (KPROMA,KLEV+1)         ; TOTAL-SKY SHORTWAVE FLUXES
+! Ajout flux LW et SW montants et descendants, et ciel clair (MPL 19.12.08)
+! PFLUX (KPROMA,2,KLEV+1)       ; LW total sky flux (1=up, 2=down)
+! PFLUC (KPROMA,2,KLEV+1)       ; LW clear sky flux (1=up, 2=down)
+! PFSDN(KPROMA,KLEV+1)          ; SW total sky flux down
+! PFSUP(KPROMA,KLEV+1)          ; SW total sky flux up
+! PFSCDN(KPROMA,KLEV+1)         ; SW clear sky flux down
+! PFSCUP(KPROMA,KLEV+1)         ; SW clear sky flux up
+
+
+!        IMPLICIT ARGUMENTS :   NONE
+!        --------------------
+
+!     METHOD.
+!     -------
+!     SEE DOCUMENTATION
+
+!     EXTERNALS.
+!     ----------
+
+!     REFERENCE.
+!     ----------
+!     ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS
+
+!     AUTHORS.
+!     --------
+!     ORIGINAL BY  B. RITTER   *ECMWF*        83-10-13
+!     REWRITING FOR IFS BY J.-J. MORCRETTE    94-11-15
+!     96-11: Ph. Dandin. Meteo-France
+!     REWRITING FOR DM  BY J.PH. PIEDELIEVRE   1998-07
+!     Duplication of RFMR to use present (cy25) ECMWF radiation scheme : Y. Bouteloup 09-2003
+!     Use of 6 aerosols & introduce NSW : F. Bouyssel 09-2004
+!     04-11-18 : 4 New arguments for AROME : Y. Seity 
+!     2005-10-10 Y. Seity : 3 optional arguments for dust optical properties 
+!     JJMorcrette 20060721 PP of clear-sky PAR and TOA incident solar radiation (ECMWF)
+!     Olivier Boucher: added LMD radiation diagnostics 2014-03
+
+!-----------------------------------------------------------------------
+
+USE PARKIND1  ,ONLY : JPIM     ,JPRB
+USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
+USE YOEAERD  , ONLY : RCAEROS
+USE YOMCST   , ONLY :         RMD      ,RMO3
+USE YOMPHY3  , ONLY : RII0
+USE YOERAD   , ONLY : NAER, RCCNLND  ,RCCNSEA  
+USE YOERDU   , ONLY : REPSCQ
+USE YOMGEM   , ONLY : NGPTOT
+USE YOERDI   , ONLY : RRAE   ,REPCLC    ,REPH2O
+USE YOMARPHY , ONLY : LRDUST 
+USE phys_output_mod, ONLY : swaero_diag
+
+!-----------------------------------------------------------------------
+
+!*       0.1   ARGUMENTS.
+!              ----------
+
+IMPLICIT NONE
+INCLUDE "clesphys.h"
+
+INTEGER(KIND=JPIM),INTENT(IN)    :: KPROMA 
+INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
+INTEGER(KIND=JPIM),INTENT(IN)    :: KST 
+INTEGER(KIND=JPIM),INTENT(IN)    :: KEND 
+INTEGER(KIND=JPIM)               :: KTDIA ! Argument NOT used
+INTEGER(KIND=JPIM),INTENT(IN)    :: KMODE
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PALBD(KPROMA,NSW) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PALBP(KPROMA,NSW) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPRS(KPROMA,KLEV+1)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPRSF(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PCCO2
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PCLFR(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PQO3(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PAER(KPROMA,KLEV,6) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PDP(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PEMIS(KPROMA) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PMU0(KPROMA) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PQ(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PQS(KPROMA,KLEV)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PQIWP(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PQLWP(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PSLM(KPROMA)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PT(KPROMA,KLEV) 
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PTS(KPROMA)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PPIZA_TOT(KPROMA,KLEV,NSW)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PCGA_TOT(KPROMA,KLEV,NSW)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAU_TOT(KPROMA,KLEV,NSW)
+!--OB
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PPIZA_NAT(KPROMA,KLEV,NSW)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PCGA_NAT(KPROMA,KLEV,NSW)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAU_NAT(KPROMA,KLEV,NSW)
+REAL(KIND=JPRB)                  :: PPIZA_ZERO(KPROMA,KLEV,NSW)
+REAL(KIND=JPRB)                  :: PCGA_ZERO(KPROMA,KLEV,NSW)
+REAL(KIND=JPRB)                  :: PTAU_ZERO(KPROMA,KLEV,NSW)
+!--fin
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_LIQ(KPROMA,KLEV)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_ICE(KPROMA,KLEV)
+!--OB
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_LIQ_PI(KPROMA,KLEV)
+REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_ICE_PI(KPROMA,KLEV)
+LOGICAL, INTENT(in)  :: ok_ade, ok_aie         ! switches whether to use aerosol direct (indirect) effects or not
+INTEGER, INTENT(in)  :: flag_aerosol           ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
+LOGICAL, INTENT(in)  :: flag_aerosol_strat     ! use stratospheric aerosols
+REAL(KIND=JPRB)   ,INTENT(out)   :: PTOPSWADAERO(KPROMA), PSOLSWADAERO(KPROMA)       ! Aerosol direct forcing at TOA and surface
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWAD0AERO(KPROMA), PSOLSWAD0AERO(KPROMA)     ! Aerosol direct forcing at TOA and surface
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWAIAERO(KPROMA), PSOLSWAIAERO(KPROMA)       ! ditto, indirect
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWCFAERO(KPROMA,3), PSOLSWCFAERO(KPROMA,3) !--do we keep this ?
+!--fin
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PEMTD(KPROMA,KLEV+1) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PEMTU(KPROMA,KLEV+1) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTRSO(KPROMA,KLEV+1) 
+REAL(KIND=JPRB)   ,INTENT(INOUT) :: PTH(KPROMA,KLEV+1)
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCTRSO(KPROMA,2) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCEMTR(KPROMA,2) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTRSOD(KPROMA) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PLWFC(KPROMA,2) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PLWFT(KPROMA,KLEV+1) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSWFC(KPROMA,2) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSWFT(KPROMA,KLEV+1) 
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSFSWDIR(KPROMA,NSW)
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSFSWDIF(KPROMA,NSW)
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSDNN(KPROMA)
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSDNV(KPROMA)
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLUX(KPROMA,2,KLEV+1) ! LW total sky flux (1=up, 2=down)
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLUC(KPROMA,2,KLEV+1) ! LW clear sky flux (1=up, 2=down)
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSDN(KPROMA,KLEV+1)   ! SW total sky flux down
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSUP(KPROMA,KLEV+1)   ! SW total sky flux up
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSCDN(KPROMA,KLEV+1)  ! SW clear sky flux down
+REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSCUP(KPROMA,KLEV+1)  ! SW clear sky flux up
+
+!     ==== COMPUTED IN RADITE ===
+!     ------------------------------------------------------------------
+!*       0.2   LOCAL ARRAYS.
+!              -------------
+REAL(KIND=JPRB) :: ZRAER  (KPROMA,6,KLEV)
+REAL(KIND=JPRB) :: ZRCLC  (KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZRMU0  (KPROMA)
+REAL(KIND=JPRB) :: ZRPR   (KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZRTI   (KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZQLWP (KPROMA,KLEV ) , ZQIWP (KPROMA,KLEV )
+
+REAL(KIND=JPRB) :: ZPQO3 (KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZQOZ (NGPTOT,KLEV)
+REAL(KIND=JPRB) :: ZQS    (KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZQ     (KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZEMTD  (KPROMA,KLEV+1)
+REAL(KIND=JPRB) :: ZEMTU  (KPROMA,KLEV+1)
+REAL(KIND=JPRB) :: ZTRSOC (KPROMA,2)
+REAL(KIND=JPRB) :: ZEMTC  (KPROMA,2)
+
+REAL(KIND=JPRB) :: ZNBAS  (KPROMA)
+REAL(KIND=JPRB) :: ZNTOP  (KPROMA)
+REAL(KIND=JPRB) :: ZQRAIN (KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZQRAINT(KPROMA,KLEV)
+REAL(KIND=JPRB) :: ZCCNL  (KPROMA)
+REAL(KIND=JPRB) :: ZCCNO  (KPROMA)
+
+!  output of radlsw
+
+REAL(KIND=JPRB) :: ZEMIT  (KPROMA)
+REAL(KIND=JPRB) :: ZFCT   (KPROMA,KLEV+1)
+REAL(KIND=JPRB) :: ZFLT   (KPROMA,KLEV+1)
+REAL(KIND=JPRB) :: ZFCS   (KPROMA,KLEV+1)
+REAL(KIND=JPRB) :: ZFLS   (KPROMA,KLEV+1)
+REAL(KIND=JPRB) :: ZFRSOD (KPROMA),ZSUDU(KPROMA)
+REAL(KIND=JPRB) :: ZPARF  (KPROMA),ZUVDF(KPROMA),ZPARCF(KPROMA),ZTINCF(KPROMA)
+
+INTEGER(KIND=JPIM) :: IBEG, IEND, JK, JL
+
+REAL(KIND=JPRB) :: ZCRAE, ZRII0, ZEMIW(KPROMA)
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+
+!---aerosol radiative diagnostics
+! Key to define the aerosol effect acting on climate
+! OB: AEROSOLFEEDBACK_ACTIVE is now a LOGICAL
+! TRUE: fluxes use natural and/or anthropogenic aerosols according to ok_ade and ok_aie, DEFAULT
+! FALSE: fluxes use no aerosols (case 1) 
+! to be used only for maintaining bit reproducibility with aerosol diagnostics activated
+LOGICAL :: AEROSOLFEEDBACK_ACTIVE = .TRUE.
+
+!OB - Fluxes including aerosol effects
+!              |        direct effect
+!ind effect    | no aerosol  NATural  TOTal
+!standard      |   5
+!natural (PI)  |               1       3     
+!total   (PD)  |               2       4   
+! so we need which case when ?
+! if flag_aerosol is on
+! ok_ade and ok_aie         = 4-2, 4-3 and 4 to proceed
+! ok_ade and not ok_aie     = 3-1 and 3 to proceed
+! not ok_ade and ok_aie     = 2-1 and 2 to proceed
+! not ok_ade and not ok_aie = 1 to proceed
+! therefore the cases have the following corresponding switches
+! 1 = not ok_ade or not ok_aie
+! 2 = ok_aie
+! 3 = ok_ade
+! 4 = ok_ade and ok_aie
+! 5 = no aerosol feedback wanted or no aerosol at all 
+! if they are called in this order then the correct call is used to proceed
+
+REAL(KIND=JPRB) ::  ZFSUP_AERO(KPROMA,KLEV+1,5)
+REAL(KIND=JPRB) ::  ZFSDN_AERO(KPROMA,KLEV+1,5)
+REAL(KIND=JPRB) ::  ZFSUP0_AERO(KPROMA,KLEV+1,5)
+REAL(KIND=JPRB) ::  ZFSDN0_AERO(KPROMA,KLEV+1,5)
+
+#include "radlsw.intfb.h"
+
+IF (LHOOK) CALL DR_HOOK('RECMWF_AERO',0,ZHOOK_HANDLE)
+IBEG=KST
+IEND=KEND
+
+!*       1.    PREPARATORY WORK
+!              ----------------
+!--OB
+!        1.0    INITIALIZATIONS
+!               --------------
+
+ZFSUP_AERO (:,:,:)=0.
+ZFSDN_AERO (:,:,:)=0.
+ZFSUP0_AERO(:,:,:)=0.
+ZFSDN0_AERO(:,:,:)=0.
+
+PTAU_ZERO(:,:,:) =1.e-15
+PPIZA_ZERO(:,:,:)=1.0
+PCGA_ZERO(:,:,:) =0.0
+
+
+!*       1.1    LOCAL CONSTANTS
+!                ---------------
+
+ZRII0=RII0
+ZCRAE=RRAE*(RRAE+2.0_JPRB)
+
+!*       2.1    FULL-LEVEL QUANTITIES
+
+ZRPR =PAPRSF
+
+DO JK=1,KLEV
+  DO JL=IBEG,IEND
+!   ZPQO3(JL,JK)=PQO3(JL,JK)*PDP(JL,JK)*RMD/RMO3
+    ZPQO3(JL,JK)=PQO3(JL,JK)*PDP(JL,JK)
+    ZRCLC(JL,JK)=MAX( 0.0_JPRB ,MIN( 1.0_JPRB ,PCLFR(JL,JK)))
+    IF (ZRCLC(JL,JK) > REPCLC) THEN
+      ZQLWP(JL,JK)=PQLWP(JL,JK)
+      ZQIWP(JL,JK)=PQIWP(JL,JK)
+    ELSE
+      ZQLWP(JL,JK)=REPH2O*ZRCLC(JL,JK)
+      ZQIWP(JL,JK)=REPH2O*ZRCLC(JL,JK)
+    ENDIF
+    ZQRAIN(JL,JK)=0.
+    ZQRAINT(JL,JK)=0.
+    ZRTI(JL,JK) =PT(JL,JK)
+    ZQS (JL,JK)=MAX(2.0_JPRB*REPH2O,PQS(JL,JK))
+    ZQ  (JL,JK)=MAX(REPH2O,MIN(PQ(JL,JK),ZQS(JL,JK)*(1.0_JPRB-REPH2O)))
+    ZEMIW(JL)=PEMIS(JL)
+  ENDDO
+ENDDO
+
+IF (NAER == 0) THEN
+  ZRAER=RCAEROS
+ELSE
+  DO JK=1,KLEV
+    DO JL=IBEG,IEND
+      ZRAER(JL,1,JK)=PAER(JL,JK,1)
+      ZRAER(JL,2,JK)=PAER(JL,JK,2)
+      ZRAER(JL,3,JK)=PAER(JL,JK,3)
+      ZRAER(JL,4,JK)=PAER(JL,JK,4)
+      ZRAER(JL,5,JK)=RCAEROS
+      ZRAER(JL,6,JK)=PAER(JL,JK,6)
+    ENDDO
+  ENDDO
+ENDIF
+
+!*       2.2    HALF-LEVEL QUANTITIES
+
+DO JK=2,KLEV
+  DO JL=IBEG,IEND
+    PTH(JL,JK)=&
+     & (PT(JL,JK-1)*PAPRSF(JL,JK-1)*(PAPRSF(JL,JK)-PAPRS(JL,JK))&
+     & +PT(JL,JK)*PAPRSF(JL,JK)*(PAPRS(JL,JK)-PAPRSF(JL,JK-1)))&
+     & *(1.0_JPRB/(PAPRS(JL,JK)*(PAPRSF(JL,JK)-PAPRSF(JL,JK-1))))  
+  ENDDO
+ENDDO
+
+!*       2.3     QUANTITIES AT BOUNDARIES
+
+DO JL=IBEG,IEND
+  PTH(JL,KLEV+1)=PTS(JL)
+  PTH(JL,1)=PT(JL,1)-PAPRSF(JL,1)*(PT(JL,1)-PTH(JL,2))&
+   & /(PAPRSF(JL,1)-PAPRS(JL,2))  
+  ZNBAS(JL)=1. 
+  ZNTOP(JL)=1.
+  ZCCNL(JL)=RCCNLND
+  ZCCNO(JL)=RCCNSEA
+ENDDO
+
+!*       3.1     SOLAR ZENITH ANGLE IS EARTH'S CURVATURE
+!                CORRECTED
+
+! CCMVAL: on impose ZRMU0=PMU0 MPL 25032010
+! 2eme essai en 3D MPL 20052010
+!DO JL=IBEG,IEND
+! ZRMU0(JL)=PMU0(JL)
+!ENDDO
+!!!!! A REVOIR MPL 20091201: enleve cette correction pour comparer a AR4
+ DO JL=IBEG,IEND
+   IF (PMU0(JL) > 1.E-10_JPRB) THEN
+     ZRMU0(JL)=RRAE/(SQRT(PMU0(JL)**2+ZCRAE)-PMU0(JL))
+   ELSE
+     ZRMU0(JL)= RRAE/SQRT(ZCRAE)
+   ENDIF   
+ ENDDO   
+
+!*         4.1     CALL TO ACTUAL RADIATION SCHEME
+!
+!----now we make multiple calls to the radiation according to which 
+!----aerosol flags are on
+
+IF (flag_aerosol .GT. 0 .OR. flag_aerosol_strat) THEN
+
+IF ( .not. ok_ade .or. .not. ok_aie ) THEN
+
+! natural aerosols for direct and indirect effect 
+! PI cloud optical properties
+! use PREF_LIQ_PI and PREF_ICE_PI
+! use NAT aerosol optical properties
+! store fluxes in index 1
+
+CALL RADLSW (&
+ & IBEG  , IEND   , KPROMA  , KLEV  , KMODE , NAER,&
+ & ZRII0 ,&
+ & ZRAER , PALBD  , PALBP   , PAPRS , ZRPR  ,&
+ & ZCCNL , ZCCNO  ,&
+ & PCCO2 , ZRCLC  , PDP     , PEMIS , ZEMIW ,PSLM    , ZRMU0 , ZPQO3,&
+ & ZQ    , ZQIWP  , ZQLWP   , ZQS   , ZQRAIN,ZQRAINT ,&
+ & PTH   , ZRTI   , PTS     , ZNBAS , ZNTOP ,&
+ & PREF_LIQ_PI, PREF_ICE_PI,&
+ & ZEMIT , ZFCT   , ZFLT    , ZFCS    , ZFLS  ,&
+ & ZFRSOD, ZSUDU  , ZUVDF   , ZPARF   , ZPARCF, ZTINCF, PSFSWDIR,&
+ & PSFSWDIF,PFSDNN, PFSDNV  ,&  
+ & LRDUST,PPIZA_NAT,PCGA_NAT,PTAU_NAT,PFLUX,PFLUC,&
+ & PFSDN , PFSUP  , PFSCDN  , PFSCUP )
+
+!* SAVE VARIABLES IN INTERIM VARIABLES A LA SW_AEROAR4
+ZFSUP0_AERO(:,:,1) = PFSCUP(:,:)
+ZFSDN0_AERO(:,:,1) = PFSCDN(:,:)
+
+ZFSUP_AERO(:,:,1) =  PFSUP(:,:)
+ZFSDN_AERO(:,:,1) =  PFSDN(:,:)
+
+ENDIF 
+
+IF (ok_aie) THEN
+
+! natural aerosols for direct indirect effect 
+! use NAT aerosol optical properties 
+! PD cloud optical properties
+! use PREF_LIQ and PREF_ICE
+! store fluxes in index 2
+
+CALL RADLSW (&
+ & IBEG  , IEND   , KPROMA  , KLEV  , KMODE , NAER,&
+ & ZRII0 ,&
+ & ZRAER , PALBD  , PALBP   , PAPRS , ZRPR  ,&
+ & ZCCNL , ZCCNO  ,&
+ & PCCO2 , ZRCLC  , PDP     , PEMIS , ZEMIW ,PSLM    , ZRMU0 , ZPQO3,&
+ & ZQ    , ZQIWP  , ZQLWP   , ZQS   , ZQRAIN,ZQRAINT ,&
+ & PTH   , ZRTI   , PTS     , ZNBAS , ZNTOP ,&
+ & PREF_LIQ, PREF_ICE,&
+ & ZEMIT , ZFCT   , ZFLT    , ZFCS    , ZFLS  ,&
+ & ZFRSOD, ZSUDU  , ZUVDF   , ZPARF   , ZPARCF, ZTINCF, PSFSWDIR,&
+ & PSFSWDIF,PFSDNN, PFSDNV  ,&  
+ & LRDUST,PPIZA_NAT,PCGA_NAT,PTAU_NAT,PFLUX,PFLUC,&
+ & PFSDN , PFSUP  , PFSCDN  , PFSCUP )
+
+!* SAVE VARIABLES IN INTERIM VARIABLES A LA SW_AEROAR4
+ZFSUP0_AERO(:,:,2) = PFSCUP(:,:)
+ZFSDN0_AERO(:,:,2) = PFSCDN(:,:)
+
+ZFSUP_AERO(:,:,2) =  PFSUP(:,:)
+ZFSDN_AERO(:,:,2) =  PFSDN(:,:)
+
+ENDIF ! ok_aie      
+
+IF (ok_ade) THEN
+
+! direct effect of total aerosol activated
+! TOT aerosols for direct effect
+! PI cloud optical properties
+! use PREF_LIQ_PI and PREF_ICE_PI
+! STORE fluxes in index 3
+ 
+CALL RADLSW (&
+ & IBEG  , IEND   , KPROMA  , KLEV  , KMODE , NAER,&
+ & ZRII0 ,&
+ & ZRAER , PALBD  , PALBP   , PAPRS , ZRPR  ,&
+ & ZCCNL , ZCCNO  ,&
+ & PCCO2 , ZRCLC  , PDP     , PEMIS , ZEMIW ,PSLM    , ZRMU0 , ZPQO3,&
+ & ZQ    , ZQIWP  , ZQLWP   , ZQS   , ZQRAIN,ZQRAINT ,&
+ & PTH   , ZRTI   , PTS     , ZNBAS , ZNTOP ,&
+ & PREF_LIQ_PI, PREF_ICE_PI,&
+ & ZEMIT , ZFCT   , ZFLT    , ZFCS    , ZFLS  ,&
+ & ZFRSOD, ZSUDU  , ZUVDF   , ZPARF   , ZPARCF, ZTINCF, PSFSWDIR,&
+ & PSFSWDIF,PFSDNN, PFSDNV  ,&  
+ & LRDUST,PPIZA_TOT,PCGA_TOT,PTAU_TOT,PFLUX,PFLUC,&
+ & PFSDN , PFSUP  , PFSCDN  , PFSCUP )
+
+!* SAVE VARIABLES IN INTERIM VARIABLES A LA SW_AEROAR4
+ZFSUP0_AERO(:,:,3) = PFSCUP(:,:)
+ZFSDN0_AERO(:,:,3) = PFSCDN(:,:)
+
+ZFSUP_AERO(:,:,3) =  PFSUP(:,:)
+ZFSDN_AERO(:,:,3) =  PFSDN(:,:)
+
+ENDIF !-end ok_ade
+
+IF (ok_ade .and. ok_aie) THEN
+
+! total aerosols for direct indirect effect 
+! use TOT aerosol optical properties 
+! PD cloud optical properties
+! use PREF_LIQ and PREF_ICE
+! store fluxes in index 4 
+
+CALL RADLSW (&
+ & IBEG  , IEND   , KPROMA  , KLEV  , KMODE , NAER,&
+ & ZRII0 ,&
+ & ZRAER , PALBD  , PALBP   , PAPRS , ZRPR  ,&
+ & ZCCNL , ZCCNO  ,&
+ & PCCO2 , ZRCLC  , PDP     , PEMIS , ZEMIW ,PSLM    , ZRMU0 , ZPQO3,&
+ & ZQ    , ZQIWP  , ZQLWP   , ZQS   , ZQRAIN,ZQRAINT ,&
+ & PTH   , ZRTI   , PTS     , ZNBAS , ZNTOP ,&
+ & PREF_LIQ, PREF_ICE,&
+ & ZEMIT , ZFCT   , ZFLT    , ZFCS    , ZFLS  ,&
+ & ZFRSOD, ZSUDU  , ZUVDF   , ZPARF   , ZPARCF, ZTINCF, PSFSWDIR,&
+ & PSFSWDIF,PFSDNN, PFSDNV  ,&  
+ & LRDUST,PPIZA_TOT,PCGA_TOT,PTAU_TOT,PFLUX,PFLUC,&
+ & PFSDN , PFSUP  , PFSCDN  , PFSCUP )
+
+!* SAVE VARIABLES IN INTERIM VARIABLES A LA SW_AEROAR4
+ZFSUP0_AERO(:,:,4) = PFSCUP(:,:)
+ZFSDN0_AERO(:,:,4) = PFSCDN(:,:)
+
+ZFSUP_AERO(:,:,4) =  PFSUP(:,:)
+ZFSDN_AERO(:,:,4) =  PFSDN(:,:)
+
+ENDIF ! ok_ade .and. ok_aie
+
+ENDIF !--if flag_aerosol GT 0 OR flag_aerosol_strat
+
+! case with no aerosols at all is also computed IF ACTIVEFEEDBACK_ACTIVE is false 
+!IF (swaero_diag .OR. .not. AEROSOLFEEDBACK_ACTIVE .OR. flag_aerosol .EQ. 0 ) THEN    
+IF (.not. AEROSOLFEEDBACK_ACTIVE .OR. flag_aerosol .EQ. 0 ) THEN    
+
+! ZERO aerosol effect
+! ZERO aerosol optical depth
+! STANDARD cloud optical properties
+! STORE fluxes in index 5
+
+CALL RADLSW (&
+ & IBEG  , IEND   , KPROMA  , KLEV  , KMODE , NAER,&
+ & ZRII0 ,&
+ & ZRAER , PALBD  , PALBP   , PAPRS , ZRPR  ,&
+ & ZCCNL , ZCCNO  ,&
+ & PCCO2 , ZRCLC  , PDP     , PEMIS , ZEMIW ,PSLM    , ZRMU0 , ZPQO3,&
+ & ZQ    , ZQIWP  , ZQLWP   , ZQS   , ZQRAIN,ZQRAINT ,&
+ & PTH   , ZRTI   , PTS     , ZNBAS , ZNTOP ,&
+!--this needs to be changed to fixed cloud optical properties
+ & PREF_LIQ_PI, PREF_ICE_PI,&
+ & ZEMIT , ZFCT   , ZFLT    , ZFCS    , ZFLS  ,&
+ & ZFRSOD, ZSUDU  , ZUVDF   , ZPARF   , ZPARCF, ZTINCF, PSFSWDIR,&
+ & PSFSWDIF,PFSDNN, PFSDNV  ,&  
+ & LRDUST,PPIZA_ZERO,PCGA_ZERO,PTAU_ZERO,PFLUX,PFLUC,&
+ & PFSDN , PFSUP  , PFSCDN  , PFSCUP )
+
+!* SAVE VARIABLES IN INTERIM VARIABLES A LA SW_AEROAR4
+ZFSUP0_AERO(:,:,5) = PFSCUP(:,:)
+ZFSDN0_AERO(:,:,5) = PFSCDN(:,:)
+
+ZFSUP_AERO(:,:,5) =  PFSUP(:,:)
+ZFSDN_AERO(:,:,5) =  PFSDN(:,:)
+
+ENDIF ! .not. AEROSOLFEEDBACK_ACTIVE
+
+!*         4.2     TRANSFORM FLUXES TO MODEL HISTORICAL VARIABLES
+
+DO JK=1,KLEV+1
+  DO JL=IBEG,IEND
+    PSWFT(JL,JK)=ZFLS(JL,JK)/(ZRII0*ZRMU0(JL))
+    PLWFT(JL,JK)=ZFLT(JL,JK)
+  ENDDO
+ENDDO
+
+ZEMTD=PLWFT
+ZEMTU=PLWFT
+
+DO JL=IBEG,IEND
+  ZTRSOC(JL, 1)=ZFCS(JL,     1)/(ZRII0*ZRMU0(JL))
+  ZTRSOC(JL, 2)=ZFCS(JL,KLEV+1)/(ZRII0*ZRMU0(JL))
+  ZEMTC (JL, 1)=ZFCT(JL,     1)
+  ZEMTC (JL, 2)=ZFCT(JL,KLEV+1)
+ENDDO
+
+!                 ------------ -- ------- -- ---- -----
+!*         5.1    STORAGE OF TRANSMISSIVITY AND EMISSIVITIES
+!*                IN KPROMA-LONG ARRAYS
+
+DO JK=1,KLEV+1
+  DO JL=IBEG,IEND
+    PEMTD(JL,JK)=ZEMTD(JL,JK)
+    PEMTU(JL,JK)=ZEMTU(JL,JK)
+    PTRSO(JL,JK)=MAX(0.0_JPRB,MIN(1.0_JPRB,PSWFT(JL,JK)))
+  ENDDO
+ENDDO
+DO JK=1,2
+  DO JL=IBEG,IEND
+    PCEMTR(JL,JK)=ZEMTC (JL,JK)
+    PCTRSO(JL,JK)=MAX( 0.0_JPRB,MIN(1.0_JPRB,ZTRSOC(JL,JK)))
+  ENDDO
+ENDDO
+DO JL=IBEG,IEND
+  PTRSOD(JL)=MAX(0.0_JPRB,MIN(1.0_JPRB,ZFRSOD(JL)/(ZRII0*ZRMU0(JL))))
+ENDDO
+
+!*         7.3   RECONSTRUCT FLUXES FOR DIAGNOSTICS
+
+DO JL=IBEG,IEND
+  IF (PMU0(JL) < 1.E-10_JPRB) ZRMU0(JL)=0.0_JPRB
+ENDDO
+DO JK=1,KLEV+1
+  DO JL=IBEG,IEND
+    PLWFT(JL,JK)=PEMTD(JL,JK)
+    PSWFT(JL,JK)=ZRMU0(JL)*ZRII0*PTRSO(JL,JK)
+  ENDDO
+ENDDO
+DO JK=1,2
+  DO JL=IBEG,IEND
+    PSWFC(JL,JK)=ZRMU0(JL)*ZRII0*PCTRSO(JL,JK)
+    PLWFC(JL,JK)=PCEMTR(JL,JK)
+  ENDDO
+ENDDO
+
+!*  8.0 DIAGNOSTICS
+!---Now we copy back the correct fields to proceed to the next timestep
+
+IF  ( AEROSOLFEEDBACK_ACTIVE .AND. (flag_aerosol .GT. 0 .OR. flag_aerosol_strat) ) THEN
+
+  IF ( ok_ade .and. ok_aie  ) THEN
+    PFSUP(:,:) =    ZFSUP_AERO(:,:,4)
+    PFSDN(:,:) =    ZFSDN_AERO(:,:,4)
+    PFSCUP(:,:) =   ZFSUP0_AERO(:,:,4)
+    PFSCDN(:,:) =   ZFSDN0_AERO(:,:,4)
+  ENDIF
+
+  IF ( ok_ade .and. (.not. ok_aie) )  THEN
+    PFSUP(:,:) =    ZFSUP_AERO(:,:,3)
+    PFSDN(:,:) =    ZFSDN_AERO(:,:,3)
+    PFSCUP(:,:) =   ZFSUP0_AERO(:,:,3)
+    PFSCDN(:,:) =   ZFSDN0_AERO(:,:,3)
+  ENDIF
+
+  IF ( (.not. ok_ade) .and. ok_aie  )  THEN
+    PFSUP(:,:) =    ZFSUP_AERO(:,:,2)
+    PFSDN(:,:) =    ZFSDN_AERO(:,:,2)
+    PFSCUP(:,:) =   ZFSUP0_AERO(:,:,2)
+    PFSCDN(:,:) =   ZFSDN0_AERO(:,:,2)
+  ENDiF
+
+  IF ((.not. ok_ade) .and. (.not. ok_aie)) THEN
+    PFSUP(:,:) =    ZFSUP_AERO(:,:,1)
+    PFSDN(:,:) =    ZFSDN_AERO(:,:,1)
+    PFSCUP(:,:) =   ZFSUP0_AERO(:,:,1)
+    PFSCDN(:,:) =   ZFSDN0_AERO(:,:,1)
+  ENDIF
+
+! The following allows to compute the forcing diagostics without
+! letting the aerosol forcing act on the meteorology
+! SEE logic above
+
+ELSE  !--not AEROSOLFEEDBACK_ACTIVE
+
+    PFSUP(:,:) =    ZFSUP_AERO(:,:,5)
+    PFSDN(:,:) =    ZFSDN_AERO(:,:,5)
+    PFSCUP(:,:) =   ZFSUP0_AERO(:,:,5)
+    PFSCDN(:,:) =   ZFSDN0_AERO(:,:,5)
+
+ENDIF
+
+!OB- HERE CHECK WITH MP IF BOTTOM AND TOP INDICES ARE OK !!!!!!!!!!!!!!!!!!
+! net anthropogenic forcing direct and 1st indirect effect diagnostics
+! requires a natural aerosol field read and used 
+! Difference of net fluxes from double call to radiation
+! Will need to be extended to LW radiation
+
+IF (flag_aerosol .GT. 0 .OR. flag_aerosol_strat) THEN
+
+IF (ok_ade.AND.ok_aie) THEN
+
+! direct anthropogenic forcing
+     PSOLSWADAERO(:)  = (ZFSDN_AERO(:,1,4)      -ZFSUP_AERO(:,1,4))      -(ZFSDN_AERO(:,1,2)      -ZFSUP_AERO(:,1,2))
+     PTOPSWADAERO(:)  = (ZFSDN_AERO(:,KLEV+1,4) -ZFSUP_AERO(:,KLEV+1,4)) -(ZFSDN_AERO(:,KLEV+1,2) -ZFSUP_AERO(:,KLEV+1,2))
+     PSOLSWAD0AERO(:) = (ZFSDN0_AERO(:,1,4)     -ZFSUP0_AERO(:,1,4))     -(ZFSDN0_AERO(:,1,2)     -ZFSUP0_AERO(:,1,2))
+     PTOPSWAD0AERO(:) = (ZFSDN0_AERO(:,KLEV+1,4)-ZFSUP0_AERO(:,KLEV+1,4))-(ZFSDN0_AERO(:,KLEV+1,2)-ZFSUP0_AERO(:,KLEV+1,2))
+
+! indirect anthropogenic forcing
+     PSOLSWAIAERO(:) = (ZFSDN_AERO(:,1,4)     -ZFSUP_AERO(:,1,4))     -(ZFSDN_AERO(:,1,3)     -ZFSUP_AERO(:,1,3))
+     PTOPSWAIAERO(:) = (ZFSDN_AERO(:,KLEV+1,4)-ZFSUP_AERO(:,KLEV+1,4))-(ZFSDN_AERO(:,KLEV+1,3)-ZFSUP_AERO(:,KLEV+1,3))
+
+! Cloud radiative forcing with natural aerosol for direct effect 
+     PSOLSWCFAERO(:,1) = (ZFSDN_AERO(:,1,2)     -ZFSUP_AERO(:,1,2))     -(ZFSDN0_AERO(:,1,2)     -ZFSUP0_AERO(:,1,2))
+     PTOPSWCFAERO(:,1) = (ZFSDN_AERO(:,KLEV+1,2)-ZFSUP_AERO(:,KLEV+1,2))-(ZFSDN0_AERO(:,KLEV+1,2)-ZFSUP0_AERO(:,KLEV+1,2))
+! Cloud radiative forcing with anthropogenic aerosol for direct effect 
+     PSOLSWCFAERO(:,2) = (ZFSDN_AERO(:,1,4)     -ZFSUP_AERO(:,1,4))     -(ZFSDN0_AERO(:,1,4)     -ZFSUP0_AERO(:,1,4))
+     PTOPSWCFAERO(:,2) = (ZFSDN_AERO(:,KLEV+1,4)-ZFSUP_AERO(:,KLEV+1,4))-(ZFSDN0_AERO(:,KLEV+1,4)-ZFSUP0_AERO(:,KLEV+1,4))
+! Cloud radiative forcing with no direct effect at all 
+     PSOLSWCFAERO(:,3) = 0.0
+     PTOPSWCFAERO(:,3) = 0.0
+
+ENDIF
+
+IF (ok_ade.AND..NOT.ok_aie) THEN
+
+! direct anthropogenic forcing
+     PSOLSWADAERO(:)  = (ZFSDN_AERO(:,1,3)      -ZFSUP_AERO(:,1,3))      -(ZFSDN_AERO(:,1,1)      -ZFSUP_AERO(:,1,1))
+     PTOPSWADAERO(:)  = (ZFSDN_AERO(:,KLEV+1,3) -ZFSUP_AERO(:,KLEV+1,3)) -(ZFSDN_AERO(:,KLEV+1,1) -ZFSUP_AERO(:,KLEV+1,1))
+     PSOLSWAD0AERO(:) = (ZFSDN0_AERO(:,1,3)     -ZFSUP0_AERO(:,1,3))     -(ZFSDN0_AERO(:,1,1)     -ZFSUP0_AERO(:,1,1))
+     PTOPSWAD0AERO(:) = (ZFSDN0_AERO(:,KLEV+1,3)-ZFSUP0_AERO(:,KLEV+1,3))-(ZFSDN0_AERO(:,KLEV+1,1)-ZFSUP0_AERO(:,KLEV+1,1))
+
+! indirect anthropogenic forcing
+     PSOLSWAIAERO(:) = 0.0
+     PTOPSWAIAERO(:) = 0.0 
+
+! Cloud radiative forcing with natural aerosol for direct effect 
+     PSOLSWCFAERO(:,1) = (ZFSDN_AERO(:,1,1)     -ZFSUP_AERO(:,1,1))     -(ZFSDN0_AERO(:,1,1)     -ZFSUP0_AERO(:,1,1))
+     PTOPSWCFAERO(:,1) = (ZFSDN_AERO(:,KLEV+1,1)-ZFSUP_AERO(:,KLEV+1,1))-(ZFSDN0_AERO(:,KLEV+1,1)-ZFSUP0_AERO(:,KLEV+1,1))
+! Cloud radiative forcing with anthropogenic aerosol for direct effect 
+     PSOLSWCFAERO(:,2) = (ZFSDN_AERO(:,1,3)     -ZFSUP_AERO(:,1,3))     -(ZFSDN0_AERO(:,1,3)     -ZFSUP0_AERO(:,1,3))
+     PTOPSWCFAERO(:,2) = (ZFSDN_AERO(:,KLEV+1,3)-ZFSUP_AERO(:,KLEV+1,3))-(ZFSDN0_AERO(:,KLEV+1,3)-ZFSUP0_AERO(:,KLEV+1,3))
+! Cloud radiative forcing with no direct effect at all 
+     PSOLSWCFAERO(:,3) = 0.0
+     PTOPSWCFAERO(:,3) = 0.0
+
+ENDIF
+
+IF (.NOT.ok_ade.AND.ok_aie) THEN
+
+! direct anthropogenic forcing
+     PSOLSWADAERO(:)  = 0.0
+     PTOPSWADAERO(:)  = 0.0
+     PSOLSWAD0AERO(:) = 0.0
+     PTOPSWAD0AERO(:) = 0.0
+
+! indirect anthropogenic forcing
+     PSOLSWAIAERO(:) = (ZFSDN_AERO(:,1,2)     -ZFSUP_AERO(:,1,2))     -(ZFSDN_AERO(:,1,1)     -ZFSUP_AERO(:,1,1))
+     PTOPSWAIAERO(:) = (ZFSDN_AERO(:,KLEV+1,2)-ZFSUP_AERO(:,KLEV+1,2))-(ZFSDN_AERO(:,KLEV+1,1)-ZFSUP_AERO(:,KLEV+1,1))
+
+! Cloud radiative forcing with natural aerosol for direct effect 
+     PSOLSWCFAERO(:,1) = (ZFSDN_AERO(:,1,2)     -ZFSUP_AERO(:,1,2))     -(ZFSDN0_AERO(:,1,2)     -ZFSUP0_AERO(:,1,2))
+     PTOPSWCFAERO(:,1) = (ZFSDN_AERO(:,KLEV+1,2)-ZFSUP_AERO(:,KLEV+1,2))-(ZFSDN0_AERO(:,KLEV+1,2)-ZFSUP0_AERO(:,KLEV+1,2))
+! Cloud radiative forcing with anthropogenic aerosol for direct effect 
+     PSOLSWCFAERO(:,2) = 0.0
+     PTOPSWCFAERO(:,2) = 0.0
+! Cloud radiative forcing with no direct effect at all 
+     PSOLSWCFAERO(:,3) = 0.0
+     PTOPSWCFAERO(:,3) = 0.0
+
+ENDIF
+
+IF (.NOT.ok_ade.AND..NOT.ok_aie) THEN
+
+! direct anthropogenic forcing
+     PSOLSWADAERO(:)  = 0.0
+     PTOPSWADAERO(:)  = 0.0
+     PSOLSWAD0AERO(:) = 0.0
+     PTOPSWAD0AERO(:) = 0.0
+
+! indirect anthropogenic forcing
+     PSOLSWAIAERO(:) = 0.0 
+     PTOPSWAIAERO(:) = 0.0
+
+! Cloud radiative forcing with natural aerosol for direct effect 
+     PSOLSWCFAERO(:,1) = (ZFSDN_AERO(:,1,1)     -ZFSUP_AERO(:,1,1))     -(ZFSDN0_AERO(:,1,1)     -ZFSUP0_AERO(:,1,1))
+     PTOPSWCFAERO(:,1) = (ZFSDN_AERO(:,KLEV+1,1)-ZFSUP_AERO(:,KLEV+1,1))-(ZFSDN0_AERO(:,KLEV+1,1)-ZFSUP0_AERO(:,KLEV+1,1))
+! Cloud radiative forcing with anthropogenic aerosol for direct effect 
+     PSOLSWCFAERO(:,2) = 0.0
+     PTOPSWCFAERO(:,2) = 0.0
+! Cloud radiative forcing with no direct effect at all 
+     PSOLSWCFAERO(:,3) = 0.0
+     PTOPSWCFAERO(:,3) = 0.0
+
+ENDIF
+
+ENDIF
+
+!IF (swaero_diag .OR. .NOT. AEROSOLFEEDBACK_ACTIVE) THEN
+IF (.NOT. AEROSOLFEEDBACK_ACTIVE) THEN
+! Cloudforcing without aerosol at all
+     PSOLSWCFAERO(:,3) = (ZFSDN_AERO(:,1,5)     -ZFSUP_AERO(:,1,5))     -(ZFSDN0_AERO(:,1,5)     -ZFSUP0_AERO(:,1,5))
+     PTOPSWCFAERO(:,3) = (ZFSDN_AERO(:,KLEV+1,5)-ZFSUP_AERO(:,KLEV+1,5))-(ZFSDN0_AERO(:,KLEV+1,5)-ZFSUP0_AERO(:,KLEV+1,5))
+ENDIF
+
+IF (LHOOK) CALL DR_HOOK('RECMWF_AERO',1,ZHOOK_HANDLE)
+END SUBROUTINE RECMWF_AERO
