Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F	(revision 1145)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F	(revision 1150)
@@ -16,12 +16,12 @@
 c Arguments:
 c
-      REAL paprs(klon,klev+1)
-      REAL pplay(klon,klev), t_seri(klon,klev)
-      REAL msulfate(klon,klev) ! masse sulfate ug SO4/m3  [ug/m^3]
-      REAL RHcl(klon,klev)     ! humidite relative ciel clair
-      REAL tau_ae(klon,klev,2) ! epaisseur optique aerosol
-      REAL piz_ae(klon,klev,2) ! single scattering albedo aerosol
-      REAL cg_ae(klon,klev,2)  ! asymmetry parameter aerosol
-      REAL ai(klon)            ! POLDER aerosol index 
+      REAL, INTENT(in) :: paprs(klon,klev+1)
+      REAL, INTENT(in) :: pplay(klon,klev), t_seri(klon,klev)
+      REAL, INTENT(in) :: msulfate(klon,klev) ! masse sulfate ug SO4/m3  [ug/m^3]
+      REAL, INTENT(in) :: RHcl(klon,klev)     ! humidite relative ciel clair
+      REAL, INTENT(out) :: tau_ae(klon,klev,2) ! epaisseur optique aerosol
+      REAL, INTENT(out) :: piz_ae(klon,klev,2) ! single scattering albedo aerosol
+      REAL, INTENT(out) :: cg_ae(klon,klev,2)  ! asymmetry parameter aerosol
+      REAL, INTENT(out) :: ai(klon)            ! POLDER aerosol index 
 c
 c Local
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90	(revision 1150)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90	(revision 1150)
@@ -0,0 +1,507 @@
+SUBROUTINE AEROPT_2BANDS( &
+     pdel, m_allaer, delt, RHcl, &
+     tau_allaer, piz_allaer, &
+     cg_allaer, fractnat_allaer, &
+     flag_aerosol, pplay, t_seri)
+
+  USE dimphy
+
+
+  !    Yves Balkanski le 12 avril 2006
+  !    Celine Deandreis
+  !    Anne Cozic Avril 2009
+  !    a partir d'une sous-routine de Johannes Quaas pour les sulfates
+  !
+  IMPLICIT NONE
+
+  INCLUDE "YOMCST.h"
+  INCLUDE "iniprint.h"
+
+  !
+  ! Input arguments:
+  !
+  REAL,                           INTENT(in)  :: pdel(KLON,KLEV)
+  REAL,                           INTENT(in)  :: delt
+  REAL, DIMENSION(klon,klev,8),   INTENT(in)  :: m_allaer
+  REAL,                           INTENT(in)  :: RHcl(KLON,KLEV)     ! humidite relative ciel clair
+  REAL, DIMENSION(klon,10),       INTENT(in)  :: fractnat_allaer
+  INTEGER,                        INTENT(in)  :: flag_aerosol
+  REAL,                           INTENT(in)  :: pplay(klon,klev)
+  REAL,                           INTENT(in)  :: t_seri(klon,klev)
+
+  !
+  ! Output arguments:
+  !
+  REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: tau_allaer ! epaisseur optique aerosol
+  REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: piz_allaer ! single scattering albedo aerosol
+  REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: cg_allaer  ! asymmetry parameter aerosol
+
+  !
+  ! Local
+  !
+  REAL, DIMENSION(klon,klev,10,2) ::  tau_ae
+  REAL, DIMENSION(klon,klev,10,2) ::  piz_ae
+  REAL, DIMENSION(klon,klev,10,2) ::  cg_ae
+  LOGICAL ::  soluble
+  INTEGER :: i, k, inu, m, mrfspecies
+  INTEGER :: spsol, spinsol
+  INTEGER :: RH_num, nbre_RH, nbsol_compaer, nbinsol_compaer
+
+  PARAMETER (nbre_RH=12)
+  PARAMETER (nbsol_compaer=5)        ! 1- Seasalt AS: 2- Sesalt CS; 3- BC soluble; 4- POM soluble; 5- SO4.
+  PARAMETER (nbinsol_compaer=3)      ! 1- Dust; 2- BC insoluble; 3- POM insoluble
+  REAL:: RH_tab(nbre_RH)
+  REAL:: RH_MAX, DELTA, rh 
+  REAL:: tau_ae2b_int(KLON,KLEV,2)   ! Intermediate computation of epaisseur optique aerosol
+  REAL:: piz_ae2b_int(KLON,KLEV,2)   ! Intermediate computation of Single scattering albedo
+  REAL:: cg_ae2b_int(KLON,KLEV,2)    ! Intermediate computation of Assymetry parameter
+  PARAMETER (RH_MAX=95.)
+  DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
+  REAL :: zrho
+  REAL :: fac
+  REAL :: zdp1(klon,klev) 
+  REAL, PARAMETER ::  gravit = 9.80616    ! m2/s
+  INTEGER, ALLOCATABLE, DIMENSION(:)  :: aerosol_name
+  INTEGER, PARAMETER :: id_SSSSM    = 1
+  INTEGER, PARAMETER :: id_CSSSM    = 2
+  INTEGER, PARAMETER :: id_ASSSM    = 3
+  INTEGER, PARAMETER :: id_ASBCM    = 4
+  INTEGER, PARAMETER :: id_ASPOMM   = 5
+  INTEGER, PARAMETER :: id_ASSO4M   = 6
+  INTEGER, PARAMETER :: id_CSSO4M   = 7
+  INTEGER, PARAMETER :: id_CIDUSTM  = 8
+  INTEGER, PARAMETER :: id_AIBCM    = 9
+  INTEGER, PARAMETER :: id_AIPOMM   = 10
+  INTEGER :: nb_aer
+  REAL, DIMENSION(klon,klev,8) :: mass_temp
+
+  !
+  ! Proprietes optiques
+  !
+  REAL:: alpha_aers_2bands(nbre_RH,2,nbsol_compaer)   !--unit m2/g SO4
+  REAL:: alpha_aeri_2bands(2,nbinsol_compaer)
+  REAL:: cg_aers_2bands(nbre_RH,2,nbsol_compaer)      !--unit 
+  REAL:: cg_aeri_2bands(2,nbinsol_compaer)
+  REAL:: piz_aers_2bands(nbre_RH,2,nbsol_compaer)     !-- unit
+  REAL:: piz_aeri_2bands(2,nbinsol_compaer)           !-- unit
+
+
+  spsol = 0
+  spinsol = 0 
+
+
+  DATA alpha_aers_2bands/  & 
+       ! seasalt soluble Coarse Soluble (CS)
+       0.5090,0.6554,0.7129,0.7767,0.8529,1.2728, & 
+       1.3820,1.5792,1.9173,2.2002,2.7173,4.1487, & 
+       0.5167,0.6613,0.7221,0.7868,0.8622,1.3027, & 
+       1.4227,1.6317,1.9887,2.2883,2.8356,4.3453, & 
+       ! seasalt soluble Accumulation Soluble (AS)
+       4.125, 4.674, 5.005, 5.434, 5.985, 10.006, & 
+       11.175,13.376,17.264,20.540,26.604, 42.349,& 
+       4.187, 3.939, 3.919, 3.937, 3.995,  5.078, & 
+       5.511, 6.434, 8.317,10.152,14.024, 26.537, &
+       ! bc soluble
+       7.675,7.675,7.675,7.675,7.675,7.675,    &
+       7.675,7.675,10.433,11.984,13.767,15.567,& 
+       4.720,4.720,4.720,4.720,4.720,4.720,    & 
+       4.720,4.720,6.081,6.793,7.567,9.344,    & 
+       ! pom soluble
+       5.503,5.503,5.503,5.503,5.588,5.957,    & 
+       6.404,7.340,8.545,10.319,13.595,20.398, & 
+       1.402,1.402,1.402,1.402,1.431,1.562,    & 
+       1.715,2.032,2.425,2.991,4.193,7.133,    & 
+       ! sulfate    
+       4.681,5.062,5.460,5.798,6.224,6.733,    & 
+       7.556,8.613,10.687,12.265,16.32,21.692, & 
+       1.107,1.239,1.381,1.490,1.635,1.8030,   &
+       2.071,2.407,3.126,3.940,5.539,7.921/ 
+
+
+  DATA alpha_aeri_2bands/  & 
+       ! dust insoluble
+       0.7661,0.7123,&
+       ! bc insoluble
+       10.360,4.437, &
+       ! pom insoluble
+       3.741,0.606/
+
+
+  DATA cg_aers_2bands/ &
+       ! seasalt Coarse soluble (CS)
+       0.727, 0.747, 0.755, 0.761, 0.770, 0.788, &
+       0.792, 0.799, 0.805, 0.809, 0.815, 0.826, &
+       0.717, 0.738, 0.745, 0.752, 0.761, 0.779, &
+       0.781, 0.786, 0.793, 0.797, 0.803, 0.813, &
+       ! Sesalt Accumulation Soluble (AS)
+       0.727, 0.741, 0.748, 0.754, 0.761, 0.782, & 
+       0.787, 0.792, 0.797, 0.799, 0.801, 0.799, &
+       0.606, 0.645, 0.658, 0.669, 0.681, 0.726, &
+       0.734, 0.746, 0.761, 0.770, 0.782, 0.798, &
+       ! bc soluble
+       .612, .612, .612, .612, .612, .612, &
+       .612, .612, .702, .734, .760, .796, &
+       .433, .433, .433, .433, .433, .433, &
+       .433, .433, .534, .575, .613, .669, &
+       ! pom soluble
+       .663, .663, .663, .663, .666, .674, &
+       .685, .702, .718, .737, .757, .777, &
+       .544, .544, .544, .544, .547, .554, &
+       .565, .583, .604, .631, .661, .698, &
+       ! sulfate    
+       .658, .669, .680, .688, .698, .707, &
+       .719, .733, .752, .760, .773, .786, &
+       .544, .555, .565, .573, .583, .593, &
+       .610, .628, .655, .666, .692, .719/
+
+
+  DATA cg_aeri_2bands/ &
+       ! dust insoluble
+       .701, .670, &
+       ! bc insoluble
+       .471, .297, &
+       ! pom insoluble
+       .568, .365/
+
+  DATA piz_aers_2bands/&
+       ! seasalt Coarse soluble (CS)
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       0.992,0.989,0.987,0.986,0.986,0.980, &
+       0.980,0.978,0.976,0.976,0.974,0.971, &
+       ! seasalt Accumulation Soluble (AS)
+       1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+       1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
+       0.970, 0.975, 0.976, 0.977, 0.978, 0.982, &
+       0.982, 0.983, 0.984, 0.984, 0.985, 0.985, & 
+       ! bc soluble
+       .445, .445, .445, .445, .445, .445, &
+       .445, .445, .461, .480, .505, .528, &
+       .362, .362, .362, .362, .362, .362, &
+       .362, .362, .381, .405, .437, .483, &
+       ! pom soluble
+       .972, .972, .972, .972, .972, .974, &
+       .976, .979, .982, .986, .989, .992, &
+       .924, .924, .924, .924, .925, .927, &
+       .932, .938, .945, .952, .961, .970, &
+       ! sulfate
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       1.000,1.000,1.000,1.000,1.000,1.000, &
+       .992, .988, .988, .987, .986, .985,  &
+       .985, .985, .984, .984, .984, .984 /
+
+
+  DATA piz_aeri_2bands/ &
+       ! dust insoluble
+       .963, .987, &
+       ! bc insoluble
+       .395, .264, &
+       ! pom insoluble
+       .966, .859/
+
+
+  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
+
+     ENDDO
+  ENDDO
+
+  IF (flag_aerosol .EQ. 1) THEN 
+     nb_aer = 1
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASSO4M
+  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 = 5
+     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
+  ELSEIF (flag_aerosol .EQ. 5) THEN 
+     nb_aer = 4
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASBCM
+     aerosol_name(2) = id_AIBCM
+     aerosol_name(3) = id_ASPOMM
+     aerosol_name(4) = id_AIPOMM
+  ELSEIF (flag_aerosol .EQ. 6) THEN 
+     nb_aer = 3
+     ALLOCATE (aerosol_name(nb_aer)) 
+     aerosol_name(1) = id_ASSO4M      
+     aerosol_name(2) = id_ASPOMM
+     aerosol_name(3) = id_AIPOMM
+  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
+
+  tau_ae(:,:,:,:)=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) 
+
+
+  DO m=1,nb_aer   ! tau is only computed for each mass
+
+     fac=1.0
+     IF (aerosol_name(m).EQ.id_SSSSM) THEN   ! for now
+        soluble=.TRUE.
+        spsol=1
+     ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN 
+        soluble=.TRUE.
+        spsol=1
+     ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN
+        soluble=.TRUE.
+        spsol=2
+     ELSEIF (aerosol_name(m).EQ.id_ASBCM) THEN
+        soluble=.TRUE.
+        spsol=3
+     ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN 
+        soluble=.TRUE.
+        spsol=4 
+     ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR.  (aerosol_name(m).EQ.id_CSSO4M)) THEN
+        soluble=.TRUE.
+        spsol=5
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+     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
+
+
+     tau_ae2b_int(:,:,:)=0.
+     piz_ae2b_int(:,:,:)=0.
+     cg_ae2b_int(:,:,:)=0.
+
+     DO k=1, KLEV
+        DO i=1, KLON
+
+           rh=MIN(RHcl(i,k)*100.,RH_MAX)
+           RH_num = INT( rh/10. + 1.)
+
+           IF (rh.GT.85.) RH_num=10
+           IF (rh.GT.90.) RH_num=11
+           DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
+
+           DO inu=1,2
+              IF (soluble) THEN
+                 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))
+
+                 tau_ae(i,k,aerosol_name(m),inu) = mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*tau_ae2b_int(i,k,inu)*fac
+
+              ELSE
+                 tau_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,5+ 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    ! boucle sur les bandes spectrale
+        ENDDO     ! Boucle sur les points géographiques (grille horizontale)
+     ENDDO     ! Boucle sur les niveaux verticaux
+  ENDDO     ! Boucle  sur les masses de traceurs
+
+
+  DO inu=1, 2
+     DO mrfspecies=1,9
+        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)
+                 
+	 	 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)+ &
+                      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,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu))
+
+
+              ELSEIF (mrfspecies .EQ. 3) THEN             ! = natural aerosol NAT
+	 	 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)+ &
+                      tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)+ &
+                      tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)+ &
+                      tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)+ &
+                      tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)+ &
+                      tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)+ &	
+                      tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)+ &
+                      tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)+ &
+                      tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)+ &
+                      tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM)
+	         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)*fractnat_allaer(i,id_ASSO4M)*piz_ae(i,k,id_ASSO4M,inu)+ &
+                      tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)*piz_ae(i,k,id_CSSO4M,inu)+ &
+                      tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)*piz_ae(i,k,id_ASBCM,inu)+ &
+                      tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)*piz_ae(i,k,id_AIBCM,inu)+ &
+                      tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)*piz_ae(i,k,id_ASPOMM,inu)+ &
+                      tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)*piz_ae(i,k,id_AIPOMM,inu)+ &	
+                      tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)*piz_ae(i,k,id_ASSSM,inu)+ &
+                      tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)*piz_ae(i,k,id_CSSSM,inu)+ &
+                      tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)*piz_ae(i,k,id_SSSSM,inu)+ &
+                      tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM)*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)*fractnat_allaer(i,id_ASSO4M)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ &
+                      tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)*piz_ae(i,k,id_CSSO4M,inu)*cg_ae(i,k,id_CSSO4M,inu)+ &
+                      tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)*piz_ae(i,k,id_ASBCM,inu)*cg_ae(i,k,id_ASBCM,inu)+ &
+                      tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu)+ &
+                      tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)*piz_ae(i,k,id_ASPOMM,inu)*cg_ae(i,k,id_ASPOMM,inu)+ &
+                      tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu)+ &
+                      tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)*piz_ae(i,k,id_ASSSM,inu)*cg_ae(i,k,id_ASSSM,inu)+ &
+                      tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu)+ &
+                      tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu)+ &
+                      tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM)*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. 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) &
+                      +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)&
+                      +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) &
+                      +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)&
+                      +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) &
+                      +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)&
+                      +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
+
+  DEALLOCATE(aerosol_name) 
+
+END SUBROUTINE AEROPT_2BANDS
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90	(revision 1150)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90	(revision 1150)
@@ -0,0 +1,515 @@
+
+
+SUBROUTINE AEROPT_5WV(&
+   pdel, m_allaer, delt, &
+   RHcl, ai, flag_aerosol, &
+   pplay, t_seri)
+
+  USE DIMPHY
+
+  !
+  !    Yves Balkanski le 12 avril 2006
+  !    Celine Deandreis
+  !    Anne Cozic  Avril 2009
+  !    a partir d'une sous-routine de Johannes Quaas pour les sulfates
+  !
+  !
+  ! 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,8), 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 
+
+  !
+  ! Local
+  !
+  INTEGER, PARAMETER :: las = 5
+  LOGICAL :: soluble
+  
+  INTEGER :: i, k, m
+  INTEGER :: spsol, spinsol, la
+  INTEGER :: RH_num
+  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 :: nbsol_compaer=5   ! 1- Seasalt AS: 2- Sesalt CS; 3- BC soluble; 4- POM soluble; 5- SO4.
+  INTEGER, PARAMETER :: nbinsol_compaer=3 ! 1- Dust; 2- BC insoluble; 3- POM insoluble
+  REAL :: zrho
+  REAL :: RH_tab(nbre_RH)
+  REAL :: DELTA, rh 
+  REAL :: tau_ae5wv_int(KLON,KLEV,las) ! Intermediate computation of epaisseur optique aerosol
+  REAL :: piz_ae5wv_int(KLON,KLEV,las) ! Intermediate single scattering albedo aerosol
+  REAL :: cg_ae5wv_int(KLON,KLEV,las)  ! Intermediate asymmetry parameter aerosol
+  REAL, PARAMETER :: RH_MAX=95.
+  DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
+  REAL :: taue670(KLON)       ! epaisseur optique aerosol absorption 550 nm
+  REAL :: taue865(KLON)       ! epaisseur optique aerosol extinction 865 nm
+  REAL :: fac
+  REAL :: zdp1(klon,klev) 
+  REAL, PARAMETER ::  gravit = 9.80616    ! m2/s
+  INTEGER, ALLOCATABLE, DIMENSION(:)  :: aerosol_name
+  INTEGER, PARAMETER :: id_SSSSM    = 1
+  INTEGER, PARAMETER :: id_CSSSM    = 2
+  INTEGER, PARAMETER :: id_ASSSM    = 3
+  INTEGER, PARAMETER :: id_ASBCM    = 4
+  INTEGER, PARAMETER :: id_ASPOMM   = 5
+  INTEGER, PARAMETER :: id_ASSO4M   = 6
+  INTEGER, PARAMETER :: id_CSSO4M   = 7
+  INTEGER, PARAMETER :: id_CIDUSTM  = 8
+  INTEGER, PARAMETER :: id_AIBCM    = 9
+  INTEGER, PARAMETER :: id_AIPOMM   = 10
+  INTEGER :: nb_aer
+  
+  REAL :: tau3d(KLON,KLEV), piz3d(KLON,KLEV), cg3d(KLON,KLEV)
+  REAL :: abs3d(KLON,KLEV)     ! epaisseur optique d'absorption
+
+  
+  REAL :: alpha_aers_5wv(nbre_RH,las,nbsol_compaer)   ! ext. coeff. Soluble comp. units *** m2/g 
+                                                      ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4.
+  REAL :: alpha_aeri_5wv(las,nbinsol_compaer)         ! ext. coeff. Insoluble comp. 1- Dust: 2- BC; 3- POM
+  REAL :: cg_aers_5wv(nbre_RH,las,nbsol_compaer)      ! Asym. param. soluble comp. 
+                                                      ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4.
+  REAL :: cg_aeri_5wv(las,nbinsol_compaer)            ! Asym. param. insoluble comp. 1- Dust: 2- BC; 3- POM
+  REAL :: piz_aers_5wv(nbre_RH,las,nbsol_compaer)   
+                                                      ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4.
+  REAL :: piz_aeri_5wv(las,nbinsol_compaer)           ! Insoluble comp. 1- Dust: 2- BC; 3- POM
+
+  REAL, ALLOCATABLE :: TAUSUM(:,:,:)     ! Aerosol optical thickness per species
+  REAL, ALLOCATABLE :: TAU(:,:,:,:)     ! Aerosol optical thickness per species
+  REAL, DIMENSION(klon,klev,8) :: mass_temp
+  
+  !
+  ! Proprietes optiques
+  !
+  REAL :: radry = 287.054                     ! dry air mass constant
+
+  !
+  ! 
+  !
+  ! 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/ &
+     ! seasalt soluble CS
+     0.50,0.90,1.05,1.21,1.40,2.41, &
+     2.66,3.11,3.88,4.52,5.69,8.84, &
+     0.51,0.92,1.07,1.23,1.42,2.45, &
+     2.70,3.16,3.94,4.58,5.76,8.94, & 
+     0.52,0.93,1.08,1.24,1.43,2.47, &
+     2.73,3.20,3.99,4.64,5.84,9.04, &
+     0.52,0.93,1.09,1.25,1.44,2.50, &
+     2.76,3.23,4.03,4.68,5.89,9.14, &
+     0.52,0.94,1.09,1.26,1.45,2.51, &
+     2.78,3.25,4.06,4.72,5.94,9.22, &
+     ! seasalt soluble AS
+     4.28, 7.17, 8.44, 9.85,11.60,22.44,  &
+     25.34,30.54,39.38,46.52,59.33,91.77, &
+     3.40, 5.67, 6.69, 7.85, 9.32,19.03,  &
+     21.78,26.88,35.87,43.40,57.33,93.43, &
+     2.48, 4.22, 5.02, 5.94, 7.11,15.29,  &
+     17.70,22.31,30.73,38.06,52.15,90.59, &
+     1.90, 3.29, 3.94, 4.69, 5.65, 12.58, &
+     14.68,18.77,26.41,33.25,46.77,85.50, &
+     1.47, 2.59, 3.12, 3.74, 4.54, 10.42, &
+     12.24,15.82,22.66,28.91,41.54,79.33, &
+     ! bc soluble
+     7.930,7.930,7.930,7.930,7.930,7.930,     &
+     7.930,7.930,10.893,12.618,14.550,16.613, &
+     7.658,7.658,7.658,7.658,7.658,7.658,     &
+     7.658,7.658,10.351,11.879,13.642,15.510, &
+     7.195,7.195,7.195,7.195,7.195,7.195,     &
+     7.195,7.195,9.551,10.847,12.381,13.994,  &
+     6.736,6.736,6.736,6.736,6.736,6.736,     &
+     6.736,6.736,8.818,9.938,11.283,12.687,   &
+     6.277,6.277,6.277,6.277,6.277,6.277,     &
+     6.277,6.277,8.123,9.094,10.275,11.501,   &
+     ! pom soluble
+     6.676,6.676,6.676,6.676,6.710,6.934,   &
+     7.141,7.569,8.034,8.529,9.456,10.511,  &
+     5.109,5.109,5.109,5.109,5.189,5.535,   &
+     5.960,6.852,8.008,9.712,12.897,19.676, &
+     3.718,3.718,3.718,3.718,3.779,4.042,   &
+     4.364,5.052,5.956,7.314,9.896,15.688,  &
+     2.849,2.849,2.849,2.849,2.897,3.107,   &
+     3.365,3.916,4.649,5.760,7.900,12.863,  &
+     2.229,2.229,2.229,2.229,2.268,2.437,   &
+     2.645,3.095,3.692,4.608,6.391,10.633,  &
+     ! Sulfate
+     5.751,6.215,6.690,7.024,7.599,8.195,      &
+     9.156,10.355,12.660,14.823,18.908,24.508, &
+     4.320,4.675,5.052,5.375,5.787,6.274,      &
+     7.066,8.083,10.088,12.003,15.697,21.133,  &
+     3.079,3.351,3.639,3.886,4.205,4.584,      &
+     5.206,6.019,7.648,9.234,12.391,17.220,    &
+     2.336,2.552,2.781,2.979,3.236,3.540,      &
+     4.046,4.711,6.056,7.388,10.093,14.313,    &
+     1.777,1.949,2.134,2.292,2.503,2.751,      &
+     3.166,3.712,4.828,5.949,8.264,11.922/
+
+  DATA alpha_aeri_5wv/ &
+     ! dust insoluble
+     0.759, 0.770, 0.775, 0.775, 0.772, &
+     !!jb bc insoluble
+     11.536,10.033, 8.422, 7.234, 6.270, &
+     ! pom insoluble
+     5.042, 3.101, 1.890, 1.294, 0.934/
+
+  DATA cg_aers_5wv/ & 
+     ! seasalt soluble (CS)
+     0.730,0.753,0.760,0.766,0.772,0.793, &
+     0.797,0.802,0.809,0.813,0.820,0.830, &
+     0.719,0.744,0.751,0.757,0.764,0.786, &
+     0.791,0.796,0.803,0.808,0.815,0.826, &
+     0.721,0.744,0.750,0.756,0.762,0.784, &
+     0.787,0.793,0.800,0.804,0.811,0.822, &
+     0.717,0.741,0.747,0.753,0.759,0.780, &
+     0.784,0.789,0.795,0.800,0.806,0.817, &
+     0.715,0.739,0.745,0.751,0.757,0.777, & 
+     0.781,0.786,0.793,0.797,0.803,0.814, &
+     ! seasalt soluble (AS)
+     0.698,0.722,0.729,0.736,0.743,0.765, &
+     0.768,0.773,0.777,0.779,0.781,0.779, &
+     0.682,0.710,0.719,0.727,0.735,0.764, &
+     0.769,0.776,0.783,0.787,0.791,0.792, &
+     0.658,0.691,0.701,0.710,0.720,0.756, &
+     0.763,0.771,0.782,0.788,0.795,0.801, &
+     0.632,0.668,0.679,0.690,0.701,0.743, &
+     0.750,0.762,0.775,0.783,0.792,0.804, &
+     0.605,0.644,0.656,0.669,0.681,0.729, &
+     0.737,0.750,0.765,0.775,0.787,0.803, &
+     ! bc soluble
+     .651, .651, .651, .651, .651, .651, &
+     .651, .651, .738, .764, .785, .800, &
+     .597, .597, .597, .597, .597, .597, &
+     .597, .597, .695, .725, .751, .770, &
+     .543, .543, .543, .543, .543, .543, &
+     .543, .543, .650, .684, .714, .736, & 
+     .504, .504, .504, .504, .504, .504, &
+     .504, .504, .614, .651, .683, .708, & 
+     .469, .469, .469, .469, .469, .469, &
+     .469, .469, .582, .620, .655, .681, &
+     ! pom soluble
+     .679, .679, .679, .679, .683, .691, &
+     .703, .720, .736, .751, .766, .784, &
+     .656, .656, .656, .656, .659, .669, &
+     .681, .699, .717, .735, .750, .779, & 
+     .623, .623, .623, .623, .627, .637, &
+     .649, .668, .688, .709, .734, .762, &
+     .592, .592, .592, .592, .595, .605, &
+     .618, .639, .660, .682, .711, .743, &
+     .561, .561, .561, .561, .565, .575, &
+     .588, .609, .632, .656, .688, .724, &
+     ! sulfate
+     .671, .684, .697, .704, .714, .723, &
+     .734, .746, .762, .771, .781, .789, &
+     .653, .666, .678, .687, .697, .707, &
+     .719, .732, .751, .762, .775, .789, &
+     .622, .635, .648, .657, .667, .678, &
+     .691, .705, .728, .741, .758, .777, &
+     .591, .604, .617, .627, .638, .650, &
+     .664, .679, .704, .719, .739, .761, &
+     .560, .574, .587, .597, .609, .621, & 
+     .637, .653, .680, .697, .719, .745/
+  !
+
+  DATA cg_aeri_5wv/&
+     ! dust insoluble
+     0.714, 0.697, 0.688, 0.683, 0.679, &
+     ! bc insoluble
+     0.511, 0.445, 0.384, 0.342, 0.307, &
+     !c pom insoluble
+     0.596, 0.536, 0.466, 0.409, 0.359/
+  !
+  DATA piz_aers_5wv/&
+     ! seasalt soluble (CS)
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     ! seasalt soluble (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, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     ! bc soluble
+     .445, .445, .445, .445, .445, .445, &
+     .445, .445, .470, .487, .508, .531, &
+     .442, .442, .442, .442, .442, .442, &
+     .442, .442, .462, .481, .506, .533, &
+     .427, .427, .427, .427, .427, .427, &
+     .427, .427, .449, .470, .497, .526, &
+     .413, .413, .413, .413, .413, .413, &
+     .413, .413, .437, .458, .486, .516, &
+     .399, .399, .399, .399, .399, .399, &
+     .399, .399, .423, .445, .473, .506, &
+     ! pom soluble
+     .975, .975, .975, .975, .975, .977, &
+     .979, .982, .984, .987, .990, .994, &
+     .972, .972, .972, .972, .973, .974, &
+     .977, .980, .983, .986, .989, .993, &
+     .963, .963, .963, .963, .964, .966, &
+     .969, .974, .977, .982, .986, .991, &
+     .955, .955, .955, .955, .955, .958, &
+     .962, .967, .972, .977, .983, .989, &
+     .944, .944, .944, .944, .944, .948, &
+     .952, .959, .962, .972, .979, .987, &
+     ! sulfate
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000, &
+     1.000,1.000,1.000,1.000,1.000,1.000/
+  !
+  DATA piz_aeri_5wv/&
+     ! dust insoluble
+     0.944, 0.970, 0.977, 0.982, 0.987, &
+     ! bc insoluble
+     0.415, 0.387, 0.355, 0.328, 0.301, &
+     ! pom insoluble
+     0.972, 0.963, 0.943, 0.923, 0.897/
+
+
+
+  DO k=1, klev
+    DO i=1, klon
+      IF (t_seri(i,k).EQ.0) stop 'stop aeropt_5wv T '
+      IF (pplay(i,k).EQ.0) stop  'stop aeropt_5wv p '
+
+      zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
+      mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9
+
+    ENDDO
+  ENDDO
+
+
+
+  IF (flag_aerosol .EQ. 1) THEN 
+      nb_aer = 1
+      ALLOCATE (aerosol_name(nb_aer)) 
+      aerosol_name(1) = id_ASSO4M
+  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 = 5
+      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
+  ELSEIF (flag_aerosol .EQ. 5) THEN 
+      nb_aer = 4
+      ALLOCATE (aerosol_name(nb_aer)) 
+      aerosol_name(1) = id_ASBCM
+      aerosol_name(2) = id_AIBCM
+      aerosol_name(3) = id_ASPOMM
+      aerosol_name(4) = id_AIPOMM
+  ELSEIF (flag_aerosol .EQ. 6) THEN 
+      nb_aer = 3
+      ALLOCATE (aerosol_name(nb_aer)) 
+      aerosol_name(1) = id_ASSO4M      
+      aerosol_name(2) = id_ASPOMM
+      aerosol_name(3) = id_AIPOMM
+  ENDIF
+      
+  ALLOCATE (tausum(klon,las,nb_aer))
+  ALLOCATE (tau(klon,klev,las,nb_aer))
+
+
+
+
+  ! 
+  ! 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
+  
+  ai(:)=0.
+  
+  tau_ae5wv_int(:,:,:)=0.
+  piz_ae5wv_int(:,:,:)=0.
+  cg_ae5wv_int(:,:,:)=0.
+  tausum(:,:,:)=0.
+  
+  tau(:,:,:,:)=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) 
+  
+  DO m=1,nb_aer   ! tau is only computed for each mass
+    
+    fac=1.0
+    IF (aerosol_name(m).EQ.id_SSSSM) THEN   ! for now
+        soluble=.TRUE.
+        spsol=1
+    ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN 
+        soluble=.TRUE.
+        spsol=1
+    ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN
+        soluble=.TRUE.
+        spsol=2
+    ELSEIF (aerosol_name(m).EQ.id_ASBCM) THEN
+        soluble=.TRUE.
+        spsol=3
+    ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN 
+        soluble=.TRUE.
+        spsol=4 
+    ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR.  (aerosol_name(m).EQ.id_CSSO4M)) THEN
+        soluble=.TRUE.
+        spsol=5
+        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
+    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
+    
+    DO la=1,las
+      tau3d(:,:)=0.
+      piz3d(:,:)=0.
+      cg3d(:,:)=0.
+      abs3d(:,:)=0.
+      
+      DO k=1, KLEV
+        DO i=1, KLON
+          
+          rh=MIN(RHcl(i,k)*100.,RH_MAX)
+          RH_num = INT( rh/10. + 1.)
+          
+          IF (rh.GT.85.) RH_num=10
+          IF (rh.GT.90.) RH_num=11
+          DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num))
+          
+          IF (soluble) THEN
+              tau_ae5wv_int(i,k,la) = &
+                 alpha_aers_5wv(RH_num,la,spsol)+DELTA* &
+                 (alpha_aers_5wv(RH_num+1,la,spsol) - & 
+                 alpha_aers_5wv(RH_num,la,spsol))
+              
+              piz_ae5wv_int(i,k,la) = &
+                 piz_aers_5wv(RH_num,la,spsol)+DELTA* &
+                 (piz_aers_5wv(RH_num+1,la,spsol) - & 
+                 piz_aers_5wv(RH_num,la,spsol))
+              
+              cg_ae5wv_int(i,k,la) = &
+                 cg_aers_5wv(RH_num,la,spsol)+DELTA* & 
+                 (cg_aers_5wv(RH_num+1,la,spsol) - & 
+                 cg_aers_5wv(RH_num,la,spsol))
+              
+              tau3d(i,k) = &
+                 mass_temp(i,k,spsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac
+
+          ELSE
+              tau_ae5wv_int(i,k,la) = alpha_aeri_5wv(la,spinsol)
+              piz_ae5wv_int(i,k,la) = piz_aeri_5wv(la,spinsol)
+              cg_ae5wv_int(i,k,la)  = cg_aeri_5wv(la,spinsol)
+
+              tau3d(i,k) = &
+                 mass_temp(i,k,5+spinsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac
+          ENDIF
+          
+          
+        ENDDO     ! Boucle sur les points géographiques (grille horizontale)
+      ENDDO     ! Boucle sur les niveaux verticaux
+      
+      tau(:,:,la,m)=tau3d(:,:)
+      
+      DO k=1, KLEV
+        DO i=1,KLON
+          tausum(i,la,m)=tausum(i,la,m)+tau3d(i,k)
+        ENDDO
+      ENDDO
+      
+    ENDDO   ! boucle sur les longueurs d'onde
+  ENDDO     ! Boucle  sur les masses de traceurs
+
+
+  taue670(:) = SUM(tausum(:,la670,:),dim=2) 
+  taue865(:) = SUM(tausum(:,la865,:),dim=2) 
+
+  DO i=1, klon
+    ai(i)=-LOG(MAX(taue670(i),0.0001)/ &
+       MAX(taue865(i),0.0001))/LOG(670./865.)
+  ENDDO
+
+  
+END SUBROUTINE AEROPT_5WV
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90	(revision 1150)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90	(revision 1150)
@@ -0,0 +1,145 @@
+SUBROUTINE aerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, &
+     pplay, paprs, t_seri, rhcl, &
+     maerosol, maerosol_pi, &
+     tau_aero, piz_aero, cg_aero )
+  
+  USE dimphy
+  IMPLICIT NONE
+
+! Input arguments
+  LOGICAL, INTENT(IN)                      :: debut
+  LOGICAL, INTENT(IN)                      :: new_aod
+  INTEGER, INTENT(IN)                      :: flag_aerosol
+  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
+
+! Output arguments
+  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: maerosol
+  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: maerosol_pi
+  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: tau_aero
+  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: piz_aero
+  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: cg_aero
+
+! Local variables
+  REAL, DIMENSION(klon)        :: aerindex ! POLDER aerosol index 
+  REAL, DIMENSION(klon,10)     :: fractnat_allaer
+  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)   :: 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)   :: pdel
+  REAL, DIMENSION(klon,klev,8) :: m_allaer
+
+  INTEGER :: k
+  
+!
+! Read aersosols from file
+!
+  maerosol(:,:) = 0. 
+  maerosol_pi(:,:) = 0. 
+
+  bcsol(:,:) = 0. 
+  bcins(:,:) = 0. 
+  pomsol(:,:) = 0. 
+  pomins(:,:) = 0. 
+  sulfate(:,:) = 0. 
+
+! Read sulfate
+  IF ( flag_aerosol .EQ. 1 .OR. &
+       flag_aerosol .EQ. 4 .OR. &
+       flag_aerosol .EQ. 6 ) THEN 
+
+     CALL readaerosol(5,rjourvrai, debut, sulfate)
+     CALL readaerosol_preind(5,rjourvrai, &
+          debut, sulfate_pi)
+
+     maerosol(:,:) = maerosol(:,:) + sulfate(:,:) 
+     maerosol_pi(:,:) = maerosol_pi(:,:) + sulfate_pi(:,:) 
+  ENDIF
+
+
+! Read bcscol and bcins
+  IF ( flag_aerosol .EQ. 2 .OR. &
+       flag_aerosol .EQ. 4 .OR. &
+       flag_aerosol .EQ. 5 ) THEN 
+
+     ! Get bc aerosol distribution 
+     CALL readaerosol(3,rjourvrai, debut,bcsol )
+     CALL readaerosol_preind(3,rjourvrai, debut, bcsol_pi)
+
+     CALL readaerosol(7,rjourvrai, debut,bcins )
+     CALL readaerosol_preind(7,rjourvrai, debut, bcins_pi)
+
+     maerosol(:,:) = maerosol(:,:) + bcsol(:,:) 
+     maerosol_pi(:,:) = maerosol_pi(:,:) + bcsol_pi(:,:) 
+  ENDIF
+
+
+! Read pomsol and pomins
+  IF ( flag_aerosol .EQ. 3 .OR. &
+       flag_aerosol .EQ. 4 .OR. &
+       flag_aerosol .EQ. 5 .OR. &
+       flag_aerosol .EQ. 6 ) THEN
+
+     CALL readaerosol(4,rjourvrai, debut,pomsol)
+     CALL readaerosol_preind(4,rjourvrai,debut,pomsol_pi )
+
+     CALL readaerosol(8,rjourvrai, debut,pomins)
+     CALL readaerosol_preind(8,rjourvrai,debut,pomins_pi )
+
+     maerosol(:,:) = maerosol(:,:) + pomsol
+     maerosol_pi(:,:) = maerosol_pi(:,:) + pomsol_pi
+  ENDIF
+
+
+! Save all aerosols in one variable
+!
+! ACo pour couplage aerosol offline 07/04/2009
+! Tableau contenant les masses pour tous les aerosols 
+! les valeurs a zero seront a remplacer par les bons
+! tableaux lorsque les routines de lectures seront
+! ajoutees. 
+  m_allaer(:,:,1) = 0.                ! SSSSM || CSSSM
+  m_allaer(:,:,2) = 0.                ! ASSSM
+  m_allaer(:,:,3) = bcsol(:,:)        ! ASBCM
+  m_allaer(:,:,4) = pomsol(:,:)       ! ASPOMM
+  m_allaer(:,:,5) = sulfate(:,:)      ! ASSO4M || CSSO4M   
+  m_allaer(:,:,6) = 0.                ! CIDUSTM
+  m_allaer(:,:,7) = bcins(:,:)        ! AIBCM
+  m_allaer(:,:,8) = pomins(:,:)       ! AIPOMM
+
+!
+! Calculate the optical properties for the aerosols
+!
+  DO k = 1, klev
+     pdel(:,k) = paprs(:,k) - paprs (:,k+1)
+  END DO
+
+  IF (.NOT. new_aod) THEN 
+     CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &
+          tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:), aerindex)
+  ELSE
+     fractnat_allaer(:,:) = 0.
+     CALL aeropt_2bands( &
+          pdel, m_allaer, pdtphys, rhcl, & 
+          tau_aero, piz_aero, cg_aero,   &
+          fractnat_allaer, flag_aerosol, &
+          pplay, t_seri) 
+
+     CALL aeropt_5wv( &
+          pdel, m_allaer, &
+          pdtphys, rhcl, aerindex, & 
+          flag_aerosol, pplay, t_seri)
+  ENDIF
+
+END SUBROUTINE aerosol_optic
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h	(revision 1145)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h	(revision 1150)
@@ -40,5 +40,5 @@
        INTEGER lev_histhf, lev_histday, lev_histmth
        CHARACTER*4 type_run
-! aer_type: pour utiliser un fichier constant dans readsulfate 
+! aer_type: pour utiliser un fichier constant dans readaerosol 
        CHARACTER*8 :: aer_type 
        LOGICAL ok_isccp, ok_regdyn
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90	(revision 1145)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90	(revision 1150)
@@ -12,5 +12,6 @@
  &                     iflag_cldcon, &
  &                     iflag_ratqs,ratqsbas,ratqshaut, &
- &		       ok_ade, ok_aie, aerosol_couple, &
+ &                     ok_ade, ok_aie, aerosol_couple, &
+ &                     flag_aerosol, new_aod, &
  &                     bl95_b0, bl95_b1,&
  &                     iflag_thermals,nsplit_thermals,tau_thermals, &
@@ -60,4 +61,6 @@
   logical              :: ok_LES
   LOGICAL              :: ok_ade, ok_aie, aerosol_couple
+  INTEGER              :: flag_aerosol
+  LOGICAL              :: new_aod
   REAL                 :: bl95_b0, bl95_b1
   real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
@@ -71,4 +74,6 @@
   logical,SAVE        :: ok_LES_omp   
   LOGICAL,SAVE        :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp
+  INTEGER, SAVE       :: flag_aerosol_omp
+  LOGICAL, SAVE       :: new_aod_omp
   REAL,SAVE           :: bl95_b0_omp, bl95_b1_omp
   REAL,SAVE           :: freq_ISCCP_omp, ecrit_ISCCP_omp
@@ -239,9 +244,34 @@
   CALL getin('aerosol_couple',aerosol_couple_omp)
 
+!
+!Config Key  = flag_aerosol
+!Config Desc = which aerosol is use for coupled model
+!Config Def  = 1
+!Config Help = Used in physiq.F
+!
+! - flag_aerosol=1 => so4 only (defaut) 
+! - flag_aerosol=2 => bc  only 
+! - flag_aerosol=3 => pom only
+! - flag_aerosol=4 => all aerosol 
+! - flag_aerosol=5 => bcpom
+! - flag_aerosol=6 => pomsulf
+
+  flag_aerosol_omp = 1
+  CALL getin('flag_aerosol',flag_aerosol_omp)
+
+! Temporary variable for testing purpose!!
+!Config Key  = new_aod
+!Config Desc = which calcul of aeropt
+!Config Def  = false
+!Config Help = Used in physiq.F
+!
+  new_aod_omp = .true.
+  CALL getin('new_aod',new_aod_omp)
+
 ! 
 !Config Key  = aer_type 
 !Config Desc = Use a constant field for the aerosols 
 !Config Def  = scenario 
-!Config Help = Used in readsulfate.F 
+!Config Help = Used in readaerosol.F90 
 ! 
   aer_type_omp = 'scenario' 
@@ -918,5 +948,5 @@
 !Config Help =
 !
-  iflag_coupl = 0
+  iflag_coupl_omp = 0
   call getin('iflag_coupl',iflag_coupl_omp)
 
@@ -927,5 +957,5 @@
 !Config Help = 
 !
-  iflag_clos = 1
+  iflag_clos_omp = 1
   call getin('iflag_clos',iflag_clos_omp)
 !
@@ -935,5 +965,5 @@
 !Config Help = 
 !
-  iflag_cvl_sigd = 0
+  iflag_cvl_sigd_omp = 0
   call getin('iflag_cvl_sigd',iflag_cvl_sigd_omp)
 
@@ -943,5 +973,5 @@
 !Config Help = 
 !
-  iflag_wake = 0
+  iflag_wake_omp = 0
   call getin('iflag_wake',iflag_wake_omp)
 
@@ -1266,4 +1296,6 @@
     ok_aie = ok_aie_omp
     aerosol_couple = aerosol_couple_omp
+    flag_aerosol=flag_aerosol_omp
+    new_aod=new_aod_omp
     aer_type = aer_type_omp
     bl95_b0 = bl95_b0_omp
@@ -1326,4 +1358,12 @@
        WRITE(numout,*)' ERROR version_ocean=',version_ocean,' not valid with slab ocean'
        CALL abort_gcm('conf_phys','version_ocean not valid',1)
+    END IF
+
+! Test sur new_aod. Ce flag permet de retrouver les resultats de l'AR4
+! il n'est utilisable que lors du couplage avec le SO4 seul 
+    IF (ok_ade .OR. ok_aie) THEN 
+       IF ( .NOT. new_aod .AND.  flag_aerosol .NE. 1) THEN
+          CALL abort_gcm('conf_phys','new_aod=.FALSE. not compatible avec flag_aerosol=1',1)
+       END IF
     END IF
 
@@ -1396,4 +1436,6 @@
   write(numout,*)' ok_aie = ',ok_aie
   write(numout,*)' aerosol_couple = ', aerosol_couple
+  write(numout,*)' flag_aerosol = ', flag_aerosol
+  write(numout,*)' new_aod = ', new_aod
   write(numout,*)' aer_type = ',aer_type
   write(numout,*)' bl95_b0 = ',bl95_b0
@@ -1409,5 +1451,5 @@
   write(numout,*)' type_run = ',type_run 
   write(numout,*)' ok_isccp = ',ok_isccp 
-  WRITE(numout,*)' solarlong0 = ', solarlong0
+  write(numout,*)' solarlong0 = ', solarlong0
   write(numout,*)' qsol0 = ', qsol0
   write(numout,*)' inertie_sol = ', inertie_sol
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h	(revision 1145)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h	(revision 1150)
@@ -832,19 +832,23 @@
 
        IF (ok_ade) THEN
-        IF (o_topswad%flag(iff)<=lev_files(iff)) THEN
-      CALL histwrite_phy(nid_files(iff),o_topswad%name,itau_w,topswad)
-        ENDIF
-        IF (o_solswad%flag(iff)<=lev_files(iff)) THEN
-      CALL histwrite_phy(nid_files(iff),o_solswad%name,itau_w,solswad)
-        ENDIF
+          IF (o_topswad%flag(iff)<=lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff),o_topswad%name,itau_w,
+     $            topswad_aero)
+          ENDIF
+          IF (o_solswad%flag(iff)<=lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff),o_solswad%name,itau_w,
+     $            solswad_aero)
+          ENDIF
        ENDIF
 
        IF (ok_aie) THEN
-        IF (o_topswai%flag(iff)<=lev_files(iff)) THEN
-      CALL histwrite_phy(nid_files(iff),o_topswai%name,itau_w,topswai)
-        ENDIF
-        IF (o_solswai%flag(iff)<=lev_files(iff)) THEN
-      CALL histwrite_phy(nid_files(iff),o_solswai%name,itau_w,solswai)
-        ENDIF
+          IF (o_topswai%flag(iff)<=lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff),o_topswai%name,itau_w,
+     $            topswai_aero)
+          ENDIF
+          IF (o_solswai%flag(iff)<=lev_files(iff)) THEN
+             CALL histwrite_phy(nid_files(iff),o_solswai%name,itau_w,
+     $            solswai_aero)
+          ENDIF
        ENDIF
 
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90	(revision 1145)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90	(revision 1150)
@@ -255,8 +255,4 @@
 !$OMP THREADPRIVATE(snow_con)
 !
-! sulfate_pi : SO4 aerosol concentration [ug/m3] (pre-industrial value)
-
-      REAL,SAVE,ALLOCATABLE :: sulfate_pi(:, :)
-!$OMP THREADPRIVATE(sulfate_pi)
       REAL,SAVE,ALLOCATABLE :: rlonPOS(:)
 !$OMP THREADPRIVATE(rlonPOS)
@@ -269,5 +265,5 @@
 ! ok_aie=T ->
 !       ok_ade=T -AIE=topswai-topswad
-!        ok_ade=F -AIE=topswai-topsw
+!       ok_ade=F -AIE=topswai-topsw
 !
 !topswad, solswad : Aerosol direct effect
@@ -277,29 +273,5 @@
       REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:)
 !$OMP THREADPRIVATE(topswai,solswai)
-      REAL,SAVE,ALLOCATABLE :: tau_ae(:,:,:), piz_ae(:,:,:)
-!$OMP THREADPRIVATE(tau_ae,piz_ae)
-      REAL,SAVE,ALLOCATABLE :: cg_ae(:,:,:)
-!$OMP THREADPRIVATE(cg_ae)
-
-! Les variables suivants uniquement pour un configuration avec INCA
-! topswad_inca, solswad_inca : Aerosol direct effect
-      REAL,SAVE,ALLOCATABLE :: topswad_inca(:), solswad_inca(:)
-!$OMP THREADPRIVATE(topswad_inca,solswad_inca)
-! topswad0_inca, solswad0_inca : Aerosol direct effect
-      REAL,SAVE,ALLOCATABLE :: topswad0_inca(:), solswad0_inca(:)
-!$OMP THREADPRIVATE(topswad0_inca,solswad0_inca)
-! topswai_inca, solswai_inca : Aerosol indirect effect
-      REAL,SAVE,ALLOCATABLE :: topswai_inca(:), solswai_inca(:)
-!$OMP THREADPRIVATE(topswai_inca,solswai_inca)
-      REAL,SAVE,ALLOCATABLE :: topsw_inca(:,:), solsw_inca(:,:)
-!$OMP THREADPRIVATE(topsw_inca,solsw_inca)
-      REAL,SAVE,ALLOCATABLE :: topsw0_inca(:,:), solsw0_inca(:,:)
-!$OMP THREADPRIVATE(topsw0_inca,solsw0_inca)
-      REAL,SAVE,ALLOCATABLE :: tau_inca(:,:,:,:)
-!$OMP THREADPRIVATE(tau_inca)
-      REAL,SAVE,ALLOCATABLE :: piz_inca(:,:,:,:)
-!$OMP THREADPRIVATE(piz_inca)
-      REAL,SAVE,ALLOCATABLE :: cg_inca(:,:,:,:)
-!$OMP THREADPRIVATE(cg_inca)
+
       REAL,SAVE,ALLOCATABLE :: ccm(:,:,:)
 !$OMP THREADPRIVATE(ccm)
@@ -402,6 +374,4 @@
       ALLOCATE(ibas_con(klon), itop_con(klon))
       ALLOCATE(rain_con(klon), snow_con(klon))
-!
-      ALLOCATE(sulfate_pi(klon, klev))
       ALLOCATE(rlonPOS(klon))
       ALLOCATE(newsst(klon))
@@ -409,19 +379,5 @@
       ALLOCATE(topswad(klon), solswad(klon))
       ALLOCATE(topswai(klon), solswai(klon))
-      ALLOCATE(tau_ae(klon,klev,2), piz_ae(klon,klev,2))
-      ALLOCATE(cg_ae(klon,klev,2))
-
-      IF (config_inca /= 'none') THEN
-         ALLOCATE(topswad_inca(klon), solswad_inca(klon))
-         ALLOCATE(topswad0_inca(klon), solswad0_inca(klon))
-         ALLOCATE(topswai_inca(klon), solswai_inca(klon))
-         ALLOCATE(topsw_inca(klon,9), solsw_inca(klon,9))
-         ALLOCATE(topsw0_inca(klon,9), solsw0_inca(klon,9))
-      END IF
-      ! Following 4 variables are needed only by INCA but must be
-      ! allocated as they exist in the phytrac argument list
-      ALLOCATE(tau_inca(klon,klev,9,2))
-      ALLOCATE(piz_inca(klon,klev,9,2))
-      ALLOCATE(cg_inca(klon,klev,9,2))
+
       ALLOCATE(ccm(klon,klev,2))
 
@@ -505,6 +461,4 @@
       deallocate(ibas_con, itop_con)
       deallocate(rain_con, snow_con)
-!
-      deallocate(sulfate_pi)
       deallocate(rlonPOS)
       deallocate(newsst)
@@ -513,17 +467,4 @@
       deallocate(topswai, solswai)
 
-      deallocate(tau_ae, piz_ae)
-      deallocate(cg_ae)
-
-      IF (config_inca /= 'none') THEN
-         deallocate(topswad_inca, solswad_inca)
-         deallocate(topswad0_inca, solswad0_inca)
-         deallocate(topswai_inca, solswai_inca)
-         deallocate(topsw_inca, solsw_inca)
-         deallocate(topsw0_inca, solsw0_inca)
-      END IF
-      deallocate(tau_inca)
-      deallocate(piz_inca)
-      deallocate(cg_inca)
       deallocate(ccm)
        
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1145)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1150)
@@ -50,5 +50,5 @@
 c
 c nlon----input-I-nombre de points horizontaux
-c nlev----input-I-nombre de couches verticales
+c nlev----input-I-nombre de couches verticales, doit etre egale a klev
 c debut---input-L-variable logique indiquant le premier passage
 c lafin---input-L-variable logique indiquant le dernier passage
@@ -736,5 +736,4 @@
       EXTERNAL phyetat0  ! lire l'etat initial de la physique
       EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
-      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
       EXTERNAL suphel    ! initialiser certaines constantes
       EXTERNAL transp    ! transport total de l'eau et de l'energie
@@ -1056,7 +1055,4 @@
       CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
       CHARACTER*40 tinst, tave, typeval
-cjq   Aerosol effects (Johannes Quaas, 27/11/2003)
-      REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
-
       REAL cldtaupi(klon,klev)  ! Cloud optical thickness for pre-industrial (pi) aerosols
 
@@ -1067,9 +1063,16 @@
 
       ! Aerosol optical properties
-
-      ! Aerosol optical properties by INCA model 
-      CHARACTER*4              ::    rfname(9) 
-      REAL aerindex(klon)       ! POLDER aerosol index
-     
+      CHARACTER*4, DIMENSION(9)      :: rfname 
+      REAL, DIMENSION(klon)          :: aerindex     ! POLDER aerosol index
+      REAL, DIMENSION(klon,klev)     :: maerosol     ! aerosol concentration [ug/m3]
+      REAL, DIMENSION(klon,klev)     :: maerosol_pi  ! aerosol concentration [ug/m3] (pre-industrial value)
+      REAL, DIMENSION(klon,klev,9,2) :: tau_aero, piz_aero, cg_aero
+      REAL, DIMENSION(klon)          :: topswad_aero, solswad_aero   ! diag
+      REAL, DIMENSION(klon)          :: topswai_aero, solswai_aero   ! diag
+      REAL, DIMENSION(klon)          :: topswad0_aero, solswad0_aero ! pas utilise, eventuellment pour diag
+      REAL, DIMENSION(klon,9)        :: topsw_aero, solsw_aero       ! pas utilise
+      REAL, DIMENSION(klon,9)        :: topsw0_aero, solsw0_aero     ! pas utilise
+
+
       ! Parameters
       LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
@@ -1080,5 +1083,9 @@
                                       ! false : lecture des aerosol dans un fichier
 c$OMP THREADPRIVATE(aerosol_couple)    
-
+      INTEGER, SAVE :: flag_aerosol 
+c$OMP THREADPRIVATE(flag_aerosol) 
+      LOGICAL, SAVE :: new_aod
+c$OMP THREADPRIVATE(new_aod) 
+   
 c
 c Declaration des constantes et des fonctions thermodynamiques
@@ -1106,15 +1113,15 @@
          write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
          write(lunout,*)
-     s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'
+     s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'
          write(lunout,*)
-     s  nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys 
+     s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys 
 
          write(lunout,*) 'papers, play, phi, u, v, t, omega'
-         do k=1,nlev
+         do k=1,klev
             write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
      s   u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
          enddo
          write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
-         do k=1,nlev
+         do k=1,klev
             write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
          enddo
@@ -1190,7 +1197,4 @@
          u10m(:,:)=0.
          v10m(:,:)=0.
-         piz_ae(:,:,:)=0.
-         tau_ae(:,:,:)=0.
-         cg_ae(:,:,:)=0.
          rain_con(:)=0.
          snow_con(:)=0.
@@ -1205,20 +1209,6 @@
          wmax_th(:)=0.
          tau_overturning_th(:)=0.
-         IF (config_inca /= 'none') THEN
-            tau_inca(:,:,:,:) = 0.
-            piz_inca(:,:,:,:) = 0.
-            cg_inca(:,:,:,:)  = 0.
-            ccm(:,:,:)        = 0.
-            topswai_inca(:)   = 0.
-            topswad_inca(:)   = 0.
-            topswad0_inca(:)  = 0.
-            topsw_inca(:,:)   = 0.
-            topsw0_inca(:,:)  = 0.
-            solswai_inca(:)   = 0.
-            solswad_inca(:)   = 0.
-            solswad0_inca(:)  = 0.
-            solsw_inca(:,:)   = 0.
-            solsw0_inca(:,:)  = 0.
-         END IF
+
+         IF (config_inca /= 'none') ccm(:,:,:) = 0.
 
          rnebcon0(:,:) = 0.0
@@ -1239,4 +1229,5 @@
      .                  iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
      .                  ok_ade, ok_aie, aerosol_couple, 
+     .                  flag_aerosol, new_aod,
      .                  bl95_b0, bl95_b1,
      .                  iflag_thermals,nsplit_thermals,tau_thermals,
@@ -2658,17 +2649,14 @@
 cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
       IF (ok_ade.OR.ok_aie) THEN
-       IF ( .NOT. aerosol_couple ) THEN
-         ! Get sulfate aerosol distribution
-         CALL readsulfate(rjourvrai, debut, sulfate)
-         CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
-
-         ! Calculate aerosol optical properties (Olivier Boucher)
-         CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
-     .        tau_ae, piz_ae, cg_ae, aerindex)
-       ENDIF
+         IF (.NOT. aerosol_couple)
+     &        CALL aerosol_optic(
+     &        debut, new_aod, flag_aerosol, rjourvrai, pdtphys,
+     &        pplay, paprs, t_seri, rhcl, 
+     &        maerosol, maerosol_pi, 
+     &        tau_aero, piz_aero, cg_aero )
       ELSE
-        tau_ae(:,:,:)=0.0
-	piz_ae(:,:,:)=0.0
-	cg_ae(:,:,:)=0.0 
+         tau_aero(:,:,:,:) = 0.
+         piz_aero(:,:,:,:) = 0.
+         cg_aero(:,:,:,:)  = 0.
       ENDIF
 
@@ -2847,7 +2835,7 @@
 
       IF (aerosol_couple) THEN 
-         sulfate(:,:) = ccm(:,:,1) 
-         sulfate_pi(:,:) = ccm(:,:,2) 
-      ENDIF
+         maerosol(:,:)    = ccm(:,:,1) 
+         maerosol_pi(:,:) = ccm(:,:,2) 
+      END IF
 
       if (ok_newmicro) then
@@ -2857,5 +2845,5 @@
      .            flwp, fiwp, flwc, fiwc,
      e            ok_aie,
-     e            sulfate, sulfate_pi,
+     e            maerosol, maerosol_pi,
      e            bl95_b0, bl95_b1,
      s            cldtaupi, re, fl)
@@ -2865,5 +2853,5 @@
      .            cldh, cldl, cldm, cldt, cldq,
      e            ok_aie,
-     e            sulfate, sulfate_pi,
+     e            maerosol, maerosol_pi,
      e            bl95_b0, bl95_b1,
      s            cldtaupi, re, fl)
@@ -2895,41 +2883,48 @@
       IF (aerosol_couple) THEN 
 #ifdef INCA
-      CALL radlwsw_inca 
-     e            (kdlon,kflev,dist, rmu0, fract, solaire,
-     e             paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
-     e             wo,
-     e             cldfra, cldemi, cldtau,
-     s             heat,heat0,cool,cool0,radsol,albpla,
-     s             topsw,toplw,solsw,sollw,
-     s             sollwdown,
-     s             topsw0,toplw0,solsw0,sollw0,
-     s             lwdn0, lwdn, lwup0, lwup, 
-     s             swdn0, swdn, swup0, swup,
-     e             ok_ade, ok_aie,
-     e             tau_inca, piz_inca, cg_inca,
-     s             topswad_inca, solswad_inca,
-     s             topswad0_inca, solswad0_inca,
-     s             topsw_inca, topsw0_inca,
-     s             solsw_inca, solsw0_inca,
-     e             cldtaupi,
-     s             topswai_inca, solswai_inca)
+         CALL radlwsw_inca 
+     e        (kdlon,kflev,dist, rmu0, fract, solaire,
+     e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
+     e        wo,
+     e        cldfra, cldemi, cldtau,
+     s        heat,heat0,cool,cool0,radsol,albpla,
+     s        topsw,toplw,solsw,sollw,
+     s        sollwdown,
+     s        topsw0,toplw0,solsw0,sollw0,
+     s        lwdn0, lwdn, lwup0, lwup, 
+     s        swdn0, swdn, swup0, swup,
+     e        ok_ade, ok_aie,
+     e        tau_aero, piz_aero, cg_aero,
+     s        topswad_aero, solswad_aero,
+     s        topswad0_aero, solswad0_aero,
+     s        topsw_aero, topsw0_aero,
+     s        solsw_aero, solsw0_aero,
+     e        cldtaupi,
+     s        topswai_aero, solswai_aero)
+            
 #endif
       ELSE
-      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
-     e            (dist, rmu0, fract, 
-     e             paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
-     e             wo,
-     e             cldfra, cldemi, cldtau,
-     s             heat,heat0,cool,cool0,radsol,albpla,
-     s             topsw,toplw,solsw,sollw,
-     s             sollwdown,
-     s             topsw0,toplw0,solsw0,sollw0,
-     s             lwdn0, lwdn, lwup0, lwup, 
-     s             swdn0, swdn, swup0, swup,
-     e             ok_ade, ok_aie, ! new for aerosol radiative effects
-     e             tau_ae, piz_ae, cg_ae, ! ="=
-     s             topswad, solswad, ! ="=
-     e             cldtaupi, ! ="=
-     s             topswai, solswai,zqsat,flwc,fiwc) ! ="=
+
+         CALL radlwsw_aero
+     e        (dist, rmu0, fract, solaire,
+     e        paprs, pplay,zxtsol,albsol1, albsol2, 
+     e        t_seri,q_seri,wo,
+     e        cldfra, cldemi, cldtau,
+     e        ok_ade, ok_aie,
+     e        tau_aero, piz_aero, cg_aero,
+     e        cldtaupi,new_aod,
+     s        heat,heat0,cool,cool0,radsol,albpla,
+     s        topsw,toplw,solsw,sollw,
+     s        sollwdown,
+     s        topsw0,toplw0,solsw0,sollw0,
+     s        lwdn0, lwdn, lwup0, lwup, 
+     s        swdn0, swdn, swup0, swup,
+     s        topswad_aero, solswad_aero,
+     s        topswai_aero, solswai_aero,
+     o        topswad0_aero, solswad0_aero,
+     o        topsw_aero, topsw0_aero,
+     o        solsw_aero, solsw0_aero)
+         
+
       ENDIF
       itaprad = 0
@@ -3159,5 +3154,5 @@
      I                   lafin,
      I                   nlon,
-     I                   nlev,
+     I                   klev,
      I                   dtime,
      I                   u,
@@ -3207,7 +3202,7 @@
      I                   aerosol_couple,
      I                   flxmass_w,
-     I                   tau_inca, 
-     I                   piz_inca, 
-     I                   cg_inca,
+     I                   tau_aero, 
+     I                   piz_aero, 
+     I                   cg_aero,
      I                   ccm,
      I                   rfname,
@@ -3218,5 +3213,5 @@
          print*,'Attention on met a 0 les thermiques pour phystoke'
 	 call phystokenc (
-     I                   nlon,nlev,pdtphys,rlon,rlat,
+     I                   nlon,klev,pdtphys,rlon,rlat,
      I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
      I                   fm_therm,entr_therm,
@@ -3415,11 +3410,11 @@
       write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
       write(lunout,*)
-     s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
+     s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
       write(lunout,*)
-     s  nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,
+     s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,
      s  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
      s  pctsrf(igout,is_sic)
       write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
-      do k=1,nlev
+      do k=1,klev
          write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
      s   d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
@@ -3427,10 +3422,10 @@
       enddo
       write(lunout,*) 'cool,heat'
-      do k=1,nlev
+      do k=1,klev
          write(lunout,*) cool(igout,k),heat(igout,k)
       enddo
 
       write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
-      do k=1,nlev
+      do k=1,klev
          write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
      s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
@@ -3439,5 +3434,5 @@
       write(lunout,*) 'd_ps ',d_ps(igout)
       write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
-      do k=1,nlev
+      do k=1,klev
          write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
      s  d_qx(igout,k,1),d_qx(igout,k,2)
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F	(revision 1145)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F	(revision 1150)
@@ -57,7 +57,7 @@
      I                    aerosol_couple,
      I                    flxmass_w,
-     I                    tau_inca, 
-     I                    piz_inca, 
-     I                    cg_inca,
+     I                    tau_aero, 
+     I                    piz_aero, 
+     I                    cg_aero,
      I                    ccm,
      I                    rfname,
@@ -138,7 +138,7 @@
       CHARACTER(len=8) :: solsym(nbtr)
       integer la
-      REAL              ::    tau_inca(klon,klev,9,2)
-      REAL              ::    piz_inca(klon,klev,9,2)
-      REAL              ::    cg_inca(klon,klev,9,2)
+      REAL              ::    tau_aero(klon,klev,9,2)
+      REAL              ::    piz_aero(klon,klev,9,2)
+      REAL              ::    cg_aero(klon,klev,9,2)
       character*4       ::    rfname(9) 
       REAL              ::    ccm(klon,klev,2) 
@@ -458,7 +458,7 @@
      $                 t_seri,       ! for chimiaq
      $                 rh,
-     $                 tau_inca,
-     $                 piz_inca,
-     $                 cg_inca,
+     $                 tau_aero,
+     $                 piz_aero,
+     $                 cg_aero,
      $                 rfname,
      $                 ccm,
@@ -497,7 +497,4 @@
      $                 sh,         !sh
      $                 rh,         !rh
-     $                 .false.,    !wrhstts
-     $                 hbuf,       !hbuf
-     $                 obuf,       !obuf
      $                 iip1,       !nx
      $                 jjp1,       !ny
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F	(revision 1145)
+++ 	(revision )
@@ -1,401 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE radlwsw(dist, rmu0, fract, 
-     .                  paprs, pplay,tsol,alb1, alb2, t,q,wo,
-     .                  cldfra, cldemi, cldtaupd,
-     .                  heat,heat0,cool,cool0,radsol,albpla,
-     .                  topsw,toplw,solsw,sollw,
-     .                  sollwdown,
-     .                  topsw0,toplw0,solsw0,sollw0,
-     .                  lwdn0, lwdn, lwup0, lwup,
-     .                  swdn0, swdn, swup0, swup,
-     .                  ok_ade, ok_aie,
-     .                  tau_ae, piz_ae, cg_ae,
-     .                  topswad, solswad,
-     .                  cldtaupi, topswai, solswai,qsat,flwc,fiwc)
-c      
-      USE dimphy
-      IMPLICIT none
-c======================================================================
-c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
-c Objet: interface entre le modele et les rayonnements
-c Arguments:
-c dist-----input-R- distance astronomique terre-soleil
-c rmu0-----input-R- cosinus de l'angle zenithal
-c fract----input-R- duree d'ensoleillement normalisee
-c co2_ppm--input-R- concentration du gaz carbonique (en ppm)
-c solaire--input-R- constante solaire (W/m**2)
-c paprs----input-R- pression a inter-couche (Pa)
-c pplay----input-R- pression au milieu de couche (Pa)
-c tsol-----input-R- temperature du sol (en K)
-c alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
-c alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge
-c t--------input-R- temperature (K)
-c q--------input-R- vapeur d'eau (en kg/kg)
-c wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505
-c cldfra---input-R- fraction nuageuse (entre 0 et 1)
-c cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
-c cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
-c ok_ade---input-L- apply the Aerosol Direct Effect or not?
-c ok_aie---input-L- apply the Aerosol Indirect Effect or not?
-c tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
-c cldtaupi-input-R- epaisseur optique des nuages dans le visible
-c                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
-c                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
-c                   it is needed for the diagnostics of the aerosol indirect radiative forcing      
-c
-c heat-----output-R- echauffement atmospherique (visible) (K/jour)
-c cool-----output-R- refroidissement dans l'IR (K/jour)
-c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
-c albpla---output-R- albedo planetaire (entre 0 et 1)
-c topsw----output-R- flux solaire net au sommet de l'atm.
-c toplw----output-R- ray. IR montant au sommet de l'atmosphere
-c solsw----output-R- flux solaire net a la surface
-c sollw----output-R- ray. IR montant a la surface
-c solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
-c topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
-c solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
-c topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
-c
-c ATTENTION: swai and swad have to be interpreted in the following manner:
-c ---------
-c ok_ade=F & ok_aie=F -both are zero
-c ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
-c                        indirect is zero
-c ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
-c                        direct is zero
-c ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
-c                        aerosol direct forcing is F_{AD} = topswai-topswad
-c
-      
-c======================================================================
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-cym#include "raddim.h"
-#include "YOETHF.h"
-c
-      real rmu0(klon), fract(klon), dist
-cIM   real co2_ppm
-cIM   real solaire
-#include "clesphys.h" 
-c
-      real paprs(klon,klev+1), pplay(klon,klev)
-      real alb1(klon), alb2(klon), tsol(klon)
-      real t(klon,klev), q(klon,klev), wo(klon,klev)
-      real cldfra(klon,klev), cldemi(klon,klev), cldtaupd(klon,klev)
-      real heat(klon,klev), cool(klon,klev)
-      real heat0(klon,klev), cool0(klon,klev)
-      real radsol(klon), topsw(klon), toplw(klon)
-      real solsw(klon), sollw(klon), albpla(klon)
-      real topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
-      real sollwdown(klon)
-cIM output 3D 
-      REAL*8 ZFSUP(KDLON,KFLEV+1)
-      REAL*8 ZFSDN(KDLON,KFLEV+1)
-      REAL*8 ZFSUP0(KDLON,KFLEV+1)
-      REAL*8 ZFSDN0(KDLON,KFLEV+1)
-c
-      REAL*8 ZFLUP(KDLON,KFLEV+1)
-      REAL*8 ZFLDN(KDLON,KFLEV+1)
-      REAL*8 ZFLUP0(KDLON,KFLEV+1)
-      REAL*8 ZFLDN0(KDLON,KFLEV+1)
-c
-      REAL*8 zx_alpha1, zx_alpha2
-c
-#include "YOMCST.h"
-c
-      INTEGER k, kk, i, j, iof, nb_gr
-      EXTERNAL LW_LMDAR4,SW_LMDAR4
-c
-cIM ctes ds clesphys.h  REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12
-      REAL*8 PSCT
-c
-      REAL*8 PALBD(kdlon,2), PALBP(kdlon,2)
-      REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
-      REAL*8 PPSOL(kdlon), PDP(kdlon,klev)
-      REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
-      REAL*8 PTAVE(kdlon,kflev)
-      REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev)
-      REAL*8 PAER(kdlon,kflev,5)
-      REAL*8 PCLDLD(kdlon,kflev)
-      REAL*8 PCLDLU(kdlon,kflev)
-      REAL*8 PCLDSW(kdlon,kflev)
-      REAL*8 PTAU(kdlon,2,kflev)
-      REAL*8 POMEGA(kdlon,2,kflev)
-      REAL*8 PCG(kdlon,2,kflev)
-c
-      REAL*8 zfract(kdlon), zrmu0(kdlon), zdist
-c
-      REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev)
-      REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev)
-      REAL*8 ztopsw(kdlon), ztoplw(kdlon)
-      REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
-cIM
-      REAL*8 zsollwdown(kdlon)
-c
-      REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
-      REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
-      REAL*8 zznormcp
-cIM output 3D : SWup, SWdn, LWup, LWdn
-      REAL swdn(klon,kflev+1),swdn0(klon,kflev+1)
-      REAL swup(klon,kflev+1),swup0(klon,kflev+1)
-      REAL lwdn(klon,kflev+1),lwdn0(klon,kflev+1)
-      REAL lwup(klon,kflev+1),lwup0(klon,kflev+1)
-      REAL qsat(klon,klev),flwc(klon,klev),fiwc(klon,klev)
-c-OB
-cjq the following quantities are needed for the aerosol radiative forcings
-
-      real topswad(klon), solswad(klon) ! output: aerosol direct forcing at TOA and surface
-      real topswai(klon), solswai(klon) ! output: aerosol indirect forcing atTOA and surface
-      real tau_ae(klon,klev,2), piz_ae(klon,klev,2), cg_ae(klon,klev,2) ! aerosol optical properties (see aeropt.F)
-      real cldtaupi(klon,klev)  ! cloud optical thickness for pre-industrial aerosol concentrations
-                                ! (i.e., with a smaller droplet concentrationand thus larger droplet radii)
-      logical ok_ade, ok_aie    ! switches whether to use aerosol direct (indirect) effects or not
-      real*8 tauae(kdlon,kflev,2) ! aer opt properties
-      real*8 pizae(kdlon,kflev,2)
-      real*8 cgae(kdlon,kflev,2)
-      REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
-      REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo
-      REAL*8 ztopswad(kdlon), zsolswad(kdlon) ! Aerosol direct forcing at TOAand surface
-      REAL*8 ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect
-cjq-end
-!rv
-      tauae(:,:,:)=0.
-      pizae(:,:,:)=0.
-      cgae(:,:,:)=0.
-!rv
-      
-c
-c-------------------------------------------
-      nb_gr = klon / kdlon
-      IF (nb_gr*kdlon .NE. klon) THEN
-         PRINT*, "kdlon mauvais:", klon, kdlon, nb_gr
-         CALL abort
-      ENDIF
-      IF (kflev .NE. klev) THEN
-          PRINT*, "kflev differe de klev, kflev, klev"
-          CALL abort
-      ENDIF
-c-------------------------------------------
-      DO k = 1, klev
-      DO i = 1, klon
-         heat(i,k)=0.
-         cool(i,k)=0.
-         heat0(i,k)=0.
-         cool0(i,k)=0.
-      ENDDO
-      ENDDO
-c
-      zdist = dist
-c
-cIM anciennes valeurs
-c     RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
-c
-cIM : on met RCO2, RCH4, RN2O, RCFC11 et RCFC12 dans clesphys.h /lecture ds conf_phys.F90
-c     RCH4 = 1.65E-06* 16.043/28.97
-c     RN2O = 306.E-09* 44.013/28.97
-c     RCFC11 = 280.E-12* 137.3686/28.97
-c     RCFC12 = 484.E-12* 120.9140/28.97
-cIM anciennes valeurs
-c     RCH4 = 1.72E-06* 16.043/28.97
-c     RN2O = 310.E-09* 44.013/28.97
-c
-c     PRINT*,'IMradlwsw : solaire, co2= ', solaire, co2_ppm
-      PSCT = solaire/zdist/zdist
-c
-      DO 99999 j = 1, nb_gr
-      iof = kdlon*(j-1)
-c
-      DO i = 1, kdlon
-         zfract(i) = fract(iof+i)
-         zrmu0(i) = rmu0(iof+i)
-         PALBD(i,1) = alb1(iof+i)
-!         PALBD(i,2) = alb1(iof+i)
-         PALBD(i,2) = alb2(iof+i)
-         PALBP(i,1) = alb1(iof+i)
-!         PALBP(i,2) = alb1(iof+i)
-         PALBP(i,2) = alb2(iof+i)
-cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96
-         PEMIS(i) = 1.0 
-         PVIEW(i) = 1.66
-         PPSOL(i) = paprs(iof+i,1)
-         zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2)) 
-     .             / (pplay(iof+i,1)-pplay(iof+i,2))
-         zx_alpha2 = 1.0 - zx_alpha1
-         PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
-         PTL(i,klev+1) = t(iof+i,klev)
-         PDT0(i) = tsol(iof+i) - PTL(i,1)
-      ENDDO
-      DO k = 2, kflev
-      DO i = 1, kdlon
-         PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
-      ENDDO
-      ENDDO
-      DO k = 1, kflev
-      DO i = 1, kdlon
-         PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
-         PTAVE(i,k) = t(iof+i,k)
-         PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
-         PQS(i,k) = PWV(i,k)
-c wo:    cm.atm (epaisseur en cm dans la situation standard)
-c POZON: kg/kg
-         POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968
-     .               /(paprs(iof+i,k)-paprs(iof+i,k+1))
-     .               *(paprs(iof+i,1)/101325.0)
-         PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
-         PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
-         PCLDSW(i,k) = cldfra(iof+i,k)
-         PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
-         PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
-         POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
-         POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
-         PCG(i,1,k) = 0.865
-         PCG(i,2,k) = 0.910
-c-OB
-cjq Introduced for aerosol indirect forcings.
-cjq The following values use the cloud optical thickness calculated from
-cjq present-day aerosol concentrations whereas the quantities without the
-cjq "A" at the end are for pre-industial (natural-only) aerosol concentrations
-cjq
-         PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
-         PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
-         POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
-         POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
-cjq-end
-      ENDDO
-      ENDDO
-c
-      DO k = 1, kflev+1
-      DO i = 1, kdlon
-         PPMB(i,k) = paprs(iof+i,k)/100.0
-      ENDDO
-      ENDDO
-c
-      DO kk = 1, 5
-      DO k = 1, kflev
-      DO i = 1, kdlon
-         PAER(i,k,kk) = 1.0E-15
-      ENDDO
-      ENDDO
-      ENDDO
-c-OB
-      DO k = 1, kflev
-      DO i = 1, kdlon
-        tauae(i,k,1)=tau_ae(iof+i,k,1)
-        pizae(i,k,1)=piz_ae(iof+i,k,1)
-        cgae(i,k,1) =cg_ae(iof+i,k,1)
-        tauae(i,k,2)=tau_ae(iof+i,k,2)
-        pizae(i,k,2)=piz_ae(iof+i,k,2)
-        cgae(i,k,2) =cg_ae(iof+i,k,2)
-      ENDDO
-      ENDDO
-c
-c===== si iflag_rrtm=0 ================================================
-cIM ctes ds clesphys.h   CALL LW(RCO2,RCH4,RN2O,RCFC11,RCFC12,
-cIM ctes ds clesphys.h   CALL SW(PSCT, RCO2, zrmu0, zfract,
-c      
-      if (iflag_rrtm.eq.0) then
-         CALL LW_LMDAR4(
-     .        PPMB, PDP,
-     .        PPSOL,PDT0,PEMIS,
-     .        PTL, PTAVE, PWV, POZON, PAER,
-     .        PCLDLD,PCLDLU,
-     .        PVIEW,
-     .        zcool, zcool0,
-     .        ztoplw,zsollw,ztoplw0,zsollw0,
-     .        zsollwdown,
-     .        ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
-         CALL SW_LMDAR4(PSCT, zrmu0, zfract,
-     S        PPMB, PDP,
-     S        PPSOL, PALBD, PALBP,
-     S        PTAVE, PWV, PQS, POZON, PAER,
-     S        PCLDSW, PTAU, POMEGA, PCG,
-     S        zheat, zheat0,
-     S        zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,
-     S        ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,
-     S        tauae, pizae, cgae, ! aerosol optical properties
-     s        PTAUA, POMEGAA,
-     s        ztopswad,zsolswad,ztopswai,zsolswai, ! diagnosed aerosol forcing
-     J        ok_ade, ok_aie) ! apply aerosol effects or not?
-      else
-c===== si iflag_rrtm=1, on passe dans SW via RECMWFL ===============
-          PRINT*, "Cette option ne fonctionne pas encore !!!"
-         CALL abort
-         endif   ! if(iflag_rrtm=0)
-
-c======================================================================
-      DO i = 1, kdlon
-         radsol(iof+i) = zsolsw(i) + zsollw(i)
-         topsw(iof+i) = ztopsw(i)
-         toplw(iof+i) = ztoplw(i)
-         solsw(iof+i) = zsolsw(i)
-         sollw(iof+i) = zsollw(i)
-         sollwdown(iof+i) = zsollwdown(i)
-cIM
-         DO k = 1, kflev+1
-         lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
-         lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
-         lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
-         lwup  ( iof+i,k)   = ZFLUP  ( i,k)
-         ENDDO
-c
-         topsw0(iof+i) = ztopsw0(i)
-         toplw0(iof+i) = ztoplw0(i)
-         solsw0(iof+i) = zsolsw0(i)
-         sollw0(iof+i) = zsollw0(i)
-         albpla(iof+i) = zalbpla(i)
-cIM
-         DO k = 1, kflev+1
-         swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
-         swdn  ( iof+i,k)   = ZFSDN  ( i,k)
-         swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
-         swup  ( iof+i,k)   = ZFSUP  ( i,k)
-         ENDDO !k=1, kflev+1
-      ENDDO
-cjq-transform the aerosol forcings, if they have
-cjq to be calculated
-      IF (ok_ade) THEN
-      DO i = 1, kdlon
-         topswad(iof+i) = ztopswad(i)
-         solswad(iof+i) = zsolswad(i)
-      ENDDO
-      ELSE
-      DO i = 1, kdlon
-         topswad(iof+i) = 0.0
-         solswad(iof+i) = 0.0
-      ENDDO
-      ENDIF
-      IF (ok_aie) THEN
-      DO i = 1, kdlon
-         topswai(iof+i) = ztopswai(i)
-         solswai(iof+i) = zsolswai(i)
-      ENDDO
-      ELSE
-      DO i = 1, kdlon
-         topswai(iof+i) = 0.0
-         solswai(iof+i) = 0.0
-      ENDDO
-      ENDIF
-cjq-end
-      DO k = 1, kflev
-c      DO i = 1, kdlon
-c         heat(iof+i,k) = zheat(i,k)
-c         cool(iof+i,k) = zcool(i,k)
-c         heat0(iof+i,k) = zheat0(i,k)
-c         cool0(iof+i,k) = zcool0(i,k)
-c      ENDDO
-      DO i = 1, kdlon
-C        scale factor to take into account the difference between
-C        dry air and watter vapour scpecific heat capacity
-         zznormcp=1.0+RVTMP2*PWV(i,k)
-         heat(iof+i,k) = zheat(i,k)/zznormcp
-         cool(iof+i,k) = zcool(i,k)/zznormcp
-         heat0(iof+i,k) = zheat0(i,k)/zznormcp
-         cool0(iof+i,k) = zcool0(i,k)/zznormcp
-      ENDDO
-      ENDDO
-c
-99999 CONTINUE
-      RETURN
-      END
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw_aero.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw_aero.F90	(revision 1150)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw_aero.F90	(revision 1150)
@@ -0,0 +1,413 @@
+SUBROUTINE radlwsw_aero( &
+   dist, rmu0, fract, solaire, &
+   paprs, pplay,tsol,albedo, alblw, &
+   t,q,wo,&
+   cldfra, cldemi, cldtaupd,&
+   ok_ade, ok_aie,&
+   tau_aero, piz_aero, cg_aero,&
+   cldtaupi, new_aod, &
+   heat,heat0,cool,cool0,radsol,albpla,&
+   topsw,toplw,solsw,sollw,&
+   sollwdown,&
+   topsw0,toplw0,solsw0,sollw0,&
+   lwdn0, lwdn, lwup0, lwup,&
+   swdn0, swdn, swup0, swup,&
+   topswad_aero, solswad_aero,&
+   topswai_aero, solswai_aero, &
+   topswad0_aero, solswad0_aero,&
+   topsw_aero, topsw0_aero,&
+   solsw_aero, solsw0_aero)
+
+
+
+  USE DIMPHY
+
+  IMPLICIT NONE
+  !======================================================================
+  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
+  ! Objet: interface entre le modele et les rayonnements
+  ! Arguments:
+  ! dist-----input-R- distance astronomique terre-soleil
+  ! rmu0-----input-R- cosinus de l'angle zenithal
+  ! fract----input-R- duree d'ensoleillement normalisee
+  ! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
+  ! solaire--input-R- constante solaire (W/m**2)
+  ! paprs----input-R- pression a inter-couche (Pa)
+  ! pplay----input-R- pression au milieu de couche (Pa)
+  ! tsol-----input-R- temperature du sol (en K)
+  ! albedo---input-R- albedo du sol (entre 0 et 1)
+  ! t--------input-R- temperature (K)
+  ! q--------input-R- vapeur d'eau (en kg/kg)
+  ! wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505
+  ! cldfra---input-R- fraction nuageuse (entre 0 et 1)
+  ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
+  ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
+  ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
+  ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
+  ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
+  ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
+  !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
+  !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
+  !                   it is needed for the diagnostics of the aerosol indirect radiative forcing      
+  !
+  ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
+  ! cool-----output-R- refroidissement dans l'IR (K/jour)
+  ! radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
+  ! albpla---output-R- albedo planetaire (entre 0 et 1)
+  ! topsw----output-R- flux solaire net au sommet de l'atm.
+  ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
+  ! solsw----output-R- flux solaire net a la surface
+  ! sollw----output-R- ray. IR montant a la surface
+  ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
+  ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
+  ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
+  ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
+  !
+  ! ATTENTION: swai and swad have to be interpreted in the following manner:
+  ! ---------
+  ! ok_ade=F & ok_aie=F -both are zero
+  ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
+  !                        indirect is zero
+  ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
+  !                        direct is zero
+  ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
+  !                        aerosol direct forcing is F_{AD} = topswai-topswad
+  !
+  
+  !======================================================================
+  
+  ! ====================================================================
+  ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
+  ! 1 = ZERO    
+  ! 2 = AER total    
+  ! 3 = NAT    
+  ! 4 = BC    
+  ! 5 = SO4    
+  ! 6 = POM    
+  ! 7 = DUST    
+  ! 8 = SS    
+  ! 9 = NO3    
+  ! 
+  ! ====================================================================
+  include "YOETHF.h"
+  include "YOMCST.h"
+
+! Input arguments
+  REAL,    INTENT(in)  :: solaire
+  REAL,    INTENT(in)  :: dist
+  REAL,    INTENT(in)  :: rmu0(KLON), fract(KLON)
+  REAL,    INTENT(in)  :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
+  REAL,    INTENT(in)  :: albedo(KLON), alblw(KLON), tsol(KLON)
+  REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV), wo(KLON,KLEV)
+  LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
+  REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
+  REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,9,2)                        ! aerosol optical properties (see aeropt.F)
+  REAL,    INTENT(in)  :: piz_aero(KLON,KLEV,9,2)                        ! aerosol optical properties (see aeropt.F)
+  REAL,    INTENT(in)  :: cg_aero(KLON,KLEV,9,2)                         ! aerosol optical properties (see aeropt.F)
+  REAL,    INTENT(in)  :: cldtaupi(KLON,KLEV)                            ! cloud optical thickness for pre-industrial aerosol concentrations
+  LOGICAL, INTENT(in)  :: new_aod                                        ! flag pour retrouver les resultats exacts de l'AR4 dans le cas ou l'on ne travaille qu'avec les sulfates
+
+! Output arguments
+  REAL,    INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
+  REAL,    INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
+  REAL,    INTENT(out) :: radsol(KLON), topsw(KLON), toplw(KLON)
+  REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON)
+  REAL,    INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
+  REAL,    INTENT(out) :: sollwdown(KLON)
+  REAL,    INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1)
+  REAL,    INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1)
+  REAL,    INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1)
+  REAL,    INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1)
+  REAL,    INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON)         ! output: aerosol direct forcing at TOA and surface
+  REAL,    INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON)         ! output: aerosol indirect forcing atTOA and surface
+  REAL, DIMENSION(klon), INTENT(out)    :: topswad0_aero 
+  REAL, DIMENSION(klon), INTENT(out)    :: solswad0_aero
+  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
+  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
+  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
+  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
+
+! Local variables
+  REAL*8 ZFSUP(KDLON,KFLEV+1)
+  REAL*8 ZFSDN(KDLON,KFLEV+1)
+  REAL*8 ZFSUP0(KDLON,KFLEV+1)
+  REAL*8 ZFSDN0(KDLON,KFLEV+1)
+  REAL*8 ZFLUP(KDLON,KFLEV+1)
+  REAL*8 ZFLDN(KDLON,KFLEV+1)
+  REAL*8 ZFLUP0(KDLON,KFLEV+1)
+  REAL*8 ZFLDN0(KDLON,KFLEV+1)
+  REAL*8 zx_alpha1, zx_alpha2
+  INTEGER k, kk, i, j, iof, nb_gr
+  REAL*8 PSCT
+  REAL*8 PALBD(kdlon,2), PALBP(kdlon,2)
+  REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
+  REAL*8 PPSOL(kdlon), PDP(kdlon,KLEV)
+  REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
+  REAL*8 PTAVE(kdlon,kflev)
+  REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev)
+  REAL*8 PAER(kdlon,kflev,5)
+  REAL*8 PCLDLD(kdlon,kflev)
+  REAL*8 PCLDLU(kdlon,kflev)
+  REAL*8 PCLDSW(kdlon,kflev)
+  REAL*8 PTAU(kdlon,2,kflev)
+  REAL*8 POMEGA(kdlon,2,kflev)
+  REAL*8 PCG(kdlon,2,kflev)
+  REAL*8 zfract(kdlon), zrmu0(kdlon), zdist
+  REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev)
+  REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev)
+  REAL*8 ztopsw(kdlon), ztoplw(kdlon)
+  REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
+  REAL*8 zsollwdown(kdlon)
+  REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
+  REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
+  REAL*8 zznormcp
+  REAL*8 tauaero(kdlon,kflev,9,2)                     ! aer opt properties
+  REAL*8 pizaero(kdlon,kflev,9,2)
+  REAL*8 cgaero(kdlon,kflev,9,2)
+  REAL*8 PTAUA(kdlon,2,kflev)                         ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
+  REAL*8 POMEGAA(kdlon,2,kflev)                       ! dito for single scatt albedo
+  REAL*8 ztopswadaero(kdlon), zsolswadaero(kdlon)     ! Aerosol direct forcing at TOAand surface
+  REAL*8 ztopswad0aero(kdlon), zsolswad0aero(kdlon)   ! Aerosol direct forcing at TOAand surface
+  REAL*8 ztopswaiaero(kdlon), zsolswaiaero(kdlon)     ! dito, indirect
+  REAL*8 ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
+  REAL*8 zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
+
+  ! initialisation
+  tauaero(:,:,:,:)=0.
+  pizaero(:,:,:,:)=0.
+  cgaero(:,:,:,:)=0.
+  
+  !
+  !-------------------------------------------
+  nb_gr = KLON / kdlon
+  IF (nb_gr*kdlon .NE. KLON) THEN
+      PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
+      CALL abort
+  ENDIF
+  IF (kflev .NE. KLEV) THEN
+      PRINT*, "kflev differe de KLEV, kflev, KLEV"
+      CALL abort
+  ENDIF
+  !-------------------------------------------
+  DO k = 1, KLEV
+    DO i = 1, KLON
+      heat(i,k)=0.
+      cool(i,k)=0.
+      heat0(i,k)=0.
+      cool0(i,k)=0.
+    ENDDO
+  ENDDO
+  !
+  zdist = dist
+  !
+  PSCT = solaire/zdist/zdist
+  DO j = 1, nb_gr
+    iof = kdlon*(j-1)
+    DO i = 1, kdlon
+      zfract(i) = fract(iof+i)
+      zrmu0(i) = rmu0(iof+i)
+      PALBD(i,1) = albedo(iof+i)
+      PALBD(i,2) = alblw(iof+i)
+      PALBP(i,1) = albedo(iof+i)
+      PALBP(i,2) = alblw(iof+i)
+      PEMIS(i) = 1.0 
+      PVIEW(i) = 1.66
+      PPSOL(i) = paprs(iof+i,1)
+      zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
+      zx_alpha2 = 1.0 - zx_alpha1
+      PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
+      PTL(i,KLEV+1) = t(iof+i,KLEV)
+      PDT0(i) = tsol(iof+i) - PTL(i,1)
+    ENDDO
+    DO k = 2, kflev
+      DO i = 1, kdlon
+        PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
+      ENDDO
+    ENDDO
+    DO k = 1, kflev
+      DO i = 1, kdlon
+        PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
+        PTAVE(i,k) = t(iof+i,k)
+        PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
+        PQS(i,k) = PWV(i,k)
+        ! wo:    cm.atm (epaisseur en cm dans la situation standard)
+        ! POZON: kg/kg
+        POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968 &
+           /(paprs(iof+i,k)-paprs(iof+i,k+1))&
+           *(paprs(iof+i,1)/101325.0)
+        PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
+        PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
+        PCLDSW(i,k) = cldfra(iof+i,k)
+        PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
+        PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
+        POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
+        POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
+        PCG(i,1,k) = 0.865
+        PCG(i,2,k) = 0.910
+        !-
+        ! Introduced for aerosol indirect forcings.
+        ! The following values use the cloud optical thickness calculated from
+        ! present-day aerosol concentrations whereas the quantities without the
+        ! "A" at the end are for pre-industial (natural-only) aerosol concentrations
+        !
+        PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
+        PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
+        POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
+        POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
+      ENDDO
+    ENDDO
+    !
+    DO k = 1, kflev+1
+      DO i = 1, kdlon
+        PPMB(i,k) = paprs(iof+i,k)/100.0
+      ENDDO
+    ENDDO
+    !
+    DO kk = 1, 5
+      DO k = 1, kflev
+        DO i = 1, kdlon
+          PAER(i,k,kk) = 1.0E-15
+        ENDDO
+      ENDDO
+    ENDDO
+    DO k = 1, kflev
+      DO i = 1, kdlon
+        tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
+        pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
+        cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
+        tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
+        pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
+        cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
+      ENDDO
+    ENDDO
+    !
+    !======================================================================
+    CALL LW_LMDAR4(&
+       PPMB, PDP,&
+       PPSOL,PDT0,PEMIS,&
+       PTL, PTAVE, PWV, POZON, PAER,&
+       PCLDLD,PCLDLU,&
+       PVIEW,&
+       zcool, zcool0,&
+       ztoplw,zsollw,ztoplw0,zsollw0,&
+       zsollwdown,&
+       ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
+
+
+    IF (.NOT. new_aod) THEN 
+       ! use old version
+        CALL SW_LMDAR4(PSCT, zrmu0, zfract,&
+           PPMB, PDP, &
+           PPSOL, PALBD, PALBP,&
+           PTAVE, PWV, PQS, POZON, PAER,&
+           PCLDSW, PTAU, POMEGA, PCG,&
+           zheat, zheat0,&
+           zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
+           ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
+           tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:),& 
+           PTAUA, POMEGAA,&
+           ztopswadaero,zsolswadaero,&
+           ztopswaiaero,zsolswaiaero,& 
+           ok_ade, ok_aie) 
+    ELSE
+
+        CALL SW_AERO(PSCT, zrmu0, zfract,&
+           PPMB, PDP,&
+           PPSOL, PALBD, PALBP,&
+           PTAVE, PWV, PQS, POZON, PAER,&
+           PCLDSW, PTAU, POMEGA, PCG,&
+           zheat, zheat0,&
+           zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
+           ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
+           tauaero, pizaero, cgaero, &
+           PTAUA, POMEGAA,&
+           ztopswadaero,zsolswadaero,&
+           ztopswad0aero,zsolswad0aero,&
+           ztopswaiaero,zsolswaiaero, & 
+           ztopsw_aero,ztopsw0_aero,&
+           zsolsw_aero,zsolsw0_aero,&
+           ok_ade, ok_aie) 
+    
+    ENDIF
+    !======================================================================
+    DO i = 1, kdlon
+      radsol(iof+i) = zsolsw(i) + zsollw(i)
+      topsw(iof+i) = ztopsw(i)
+      toplw(iof+i) = ztoplw(i)
+      solsw(iof+i) = zsolsw(i)
+      sollw(iof+i) = zsollw(i)
+      sollwdown(iof+i) = zsollwdown(i)
+      DO k = 1, kflev+1
+        lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
+        lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
+        lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
+        lwup  ( iof+i,k)   = ZFLUP  ( i,k)
+      ENDDO
+      topsw0(iof+i) = ztopsw0(i)
+      toplw0(iof+i) = ztoplw0(i)
+      solsw0(iof+i) = zsolsw0(i)
+      sollw0(iof+i) = zsollw0(i)
+      albpla(iof+i) = zalbpla(i)
+
+      DO k = 1, kflev+1
+        swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
+        swdn  ( iof+i,k)   = ZFSDN  ( i,k)
+        swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
+        swup  ( iof+i,k)   = ZFSUP  ( i,k)
+      ENDDO
+    ENDDO
+    !-transform the aerosol forcings, if they have
+    ! to be calculated
+    IF (ok_ade) THEN
+        DO i = 1, kdlon
+          topswad_aero(iof+i) = ztopswadaero(i)
+          topswad0_aero(iof+i) = ztopswad0aero(i)
+          solswad_aero(iof+i) = zsolswadaero(i)
+          solswad0_aero(iof+i) = zsolswad0aero(i)
+          topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
+          topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
+          solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
+          solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
+          
+        ENDDO
+    ELSE
+        DO i = 1, kdlon
+          topswad_aero(iof+i) = 0.0
+          solswad_aero(iof+i) = 0.0
+          topswad0_aero(iof+i) = 0.0
+          solswad0_aero(iof+i) = 0.0
+          topsw_aero(iof+i,:) = 0.
+          topsw0_aero(iof+i,:) =0.
+          solsw_aero(iof+i,:) = 0.
+          solsw0_aero(iof+i,:) = 0.
+        ENDDO
+    ENDIF
+    IF (ok_aie) THEN
+        DO i = 1, kdlon
+          topswai_aero(iof+i) = ztopswaiaero(i)
+          solswai_aero(iof+i) = zsolswaiaero(i)
+        ENDDO
+    ELSE
+        DO i = 1, kdlon
+          topswai_aero(iof+i) = 0.0
+          solswai_aero(iof+i) = 0.0
+        ENDDO
+    ENDIF
+    DO k = 1, kflev
+      DO i = 1, kdlon
+        !        scale factor to take into account the difference between
+        !        dry air and watter vapour scpecifi! heat capacity
+        zznormcp=1.0+RVTMP2*PWV(i,k)
+        heat(iof+i,k) = zheat(i,k)/zznormcp
+        cool(iof+i,k) = zcool(i,k)/zznormcp
+        heat0(iof+i,k) = zheat0(i,k)/zznormcp
+        cool0(iof+i,k) = zcool0(i,k)/zznormcp
+      ENDDO
+    ENDDO
+    !
+  ENDDO
+
+
+ENDSUBROUTINE radlwsw_aero
+
+
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90	(revision 1150)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90	(revision 1150)
@@ -0,0 +1,701 @@
+!
+! $Header$
+!
+SUBROUTINE readaerosol (id_aero, r_day, first, massvar_p)
+  
+  USE dimphy, ONLY : klev
+  USE mod_grid_phy_lmdz, klon=>klon_glo
+  USE mod_phys_lmdz_para
+  
+  IMPLICIT NONE
+      
+! Content: 
+! --------
+! This routine reads in monthly mean values of massvar aerosols and 
+! returns a linearly interpolated dayly-mean field.      
+! 
+!
+! Author:
+! -------
+! Johannes Quaas (quaas@lmd.jussieu.fr) 
+! Celine Deandreis & Anne Cozic
+! 19/02/09
+!
+  
+!      
+! Include-files:
+! --------------     
+  INCLUDE "YOMCST.h"
+  INCLUDE "chem.h"      
+  INCLUDE "dimensions.h"      
+  INCLUDE "temps.h"      
+  INCLUDE "clesphys.h"
+  INCLUDE "iniprint.h"
+!
+! Input:
+! ------
+  INTEGER, INTENT(in) :: id_aero
+  REAL, INTENT(in)    :: r_day        ! Day of integration
+  LOGICAL, INTENT(in) :: first        ! First timestep 
+                                      ! (and therefore initialization necessary)
+!      
+! Output:      
+! -------     
+  REAL, INTENT(out) ::   massvar_p(klon_omp,klev) ! Mass of massvar (monthly mean data,from file) [ug AIBCM/m3]
+
+!      
+! Local Variables:
+! ----------------      
+  INTEGER                         ::  i, ig, k, it
+  INTEGER                         :: j, iday, iyr, iyr1, iyr2
+  INTEGER                         :: im, day1, day2, im2
+  INTEGER, PARAMETER              :: ny=jjm+1
+  CHARACTER(len=4)                :: cyear
+  CHARACTER(len=7),DIMENSION(8)   :: name_aero
+  REAL                            :: var_1(iim, jjm+1, klev, 12)
+  REAL, DIMENSION(klon,klev)      ::  massvar
+  REAL, DIMENSION(iim, jjm+1, klev, 12)        :: var_2  ! The massvar distributions
+  REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE  :: var ! VAR in right dimension
+  REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE     :: var_out
+!$OMP THREADPRIVATE(var,var_out)
+
+  LOGICAL            :: lnewday 
+  LOGICAL, PARAMETER :: lonlyone=.FALSE.
+  LOGICAL,SAVE       :: first2=.TRUE.
+!$OMP THREADPRIVATE(first2)
+
+! Variable pour definir a partir d'un numero d'aerosol, son nom
+  name_aero(1) = "SSSSM  "
+  name_aero(2) = "ASSSM  "
+  name_aero(3) = "ASBCM  "
+  name_aero(4) = "ASPOMM "
+  name_aero(5) = "SO4    "
+  name_aero(6) = "CIDUSTM"
+  name_aero(7) = "AIBCM  "
+  name_aero(8) = "AIPOMM "
+
+!$OMP MASTER
+  IF (first2) THEN
+      ALLOCATE( var(klon, klev, 12,8) )
+      ALLOCATE( var_out(klon, klev,8))
+      first2=.FALSE.
+
+      IF (aer_type /= 'actuel  ' .AND. aer_type /= 'preind  ' .AND.   &
+           aer_type /= 'scenario') THEN
+         WRITE(lunout,*)' *** Warning ***'
+         WRITE(lunout,*)'Option aer_type non prevu pour les aerosols = ',       &
+              aer_type
+         CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
+      ENDIF
+   ENDIF
+
+  IF (is_mpi_root) THEN
+      
+      iday = INT(r_day) 
+      ! Get the year of the run
+      iyr  = iday/360
+      ! Get the day of the actual year:
+      iday = iday-iyr*360
+      ! 0.02 is about 0.5/24, namly less than half an hour
+      lnewday = (r_day-FLOAT(iday).LT.0.02)
+      
+      !     ---------------------------------------------
+      !     All has to be done only, if a new day begins!       
+      !     ---------------------------------------------
+         
+      IF (lnewday.OR.first) THEN
+          
+          im = iday/30 +1     ! the actual month
+          ! annee_ref is the initial year (defined in temps.h)
+          iyr = iyr + annee_ref
+          
+          ! Do I have to read new data? (Is this the first day of a year?)
+          IF (first.OR.iday.EQ.1.) THEN 
+              ! Initialize values
+              DO it=1,12
+                DO k=1,klev
+                  DO i=1,klon
+                    var(i,k,it,id_aero)=0.
+                  ENDDO
+                ENDDO
+              ENDDO
+              
+
+              IF (aer_type == 'actuel  ') then
+
+                 cyear='1980'
+                 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
+
+              ELSE IF (aer_type == 'preind  ') THEN
+
+                 cyear='.nat'
+                 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
+
+              ELSE ! aer_type=scenario
+
+                 IF (iyr .LT. 1850) THEN
+                    cyear='.nat'
+                    WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
+                    CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
+
+                 ELSE IF (iyr .GE. 2100) THEN
+                    cyear='2100'
+                    WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
+                    CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
+                 ELSE
+                    ! Read in data:
+                    ! a) from actual 10-yr-period
+                    IF (iyr.LT.1900) THEN
+                       iyr1 = 1850
+                       iyr2 = 1900
+                    ELSE IF (iyr.GE.1900.AND.iyr.LT.1920) THEN
+                       iyr1 = 1900
+                       iyr2 = 1920
+                    ELSE 
+                       iyr1 = INT(iyr/10)*10
+                       iyr2 = INT(1+iyr/10)*10
+                    ENDIF
+                    
+                    WRITE(cyear,'(I4)') iyr1
+                    WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
+                    
+                    CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
+                    
+                 ENDIF
+                 ! If to read two decades:
+                 IF (.NOT.lonlyone) THEN
+                    
+                    ! b) from the next following one
+                    WRITE(cyear,'(I4)') iyr2
+                    WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
+                    
+                    CALL get_aero_fromfile(cyear, var_2, name_aero(id_aero))
+                    
+                    ! Interpolate linarily to the actual year:
+                    DO it=1,12
+                       DO k=1,klev
+                          DO j=1,jjm
+                             DO i=1,iim
+                                var_1(i,j,k,it) = &
+                                     var_1(i,j,k,it) - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) * &
+                                     (var_1(i,j,k,it) - var_2(i,j,k,it))
+                             ENDDO
+                          ENDDO
+                       ENDDO
+                    ENDDO
+                    
+                 ENDIF     !lonlyone
+              ENDIF ! aer_type
+               
+              ! Transform the horizontal 2D-field into the physics-field
+              ! (Also the levels and the latitudes have to be inversed)
+              
+              DO it=1,12
+                DO k=1, klev         
+                  ! a) at the poles, use the zonal mean:
+                  DO i=1,iim
+                    ! North pole
+                    var(1,k,it,id_aero) = &
+                       var(1,k,it,id_aero)+ var_1(i,jjm+1,klev+1-k,it)
+                    ! South pole
+                    var(klon,k,it,id_aero)= &
+                       var(klon,k,it,id_aero)+ var_1(i,1,klev+1-k,it)
+                  ENDDO
+                  var(1,k,it,id_aero)   = var(1,k,it,id_aero)/FLOAT(iim)
+                  var(klon,k,it,id_aero)= var(klon,k,it,id_aero)/FLOAT(iim)
+                   
+                  ! b) the values between the poles:
+                  ig=1
+                  DO j=2,jjm
+                    DO i=1,iim
+                      ig=ig+1
+
+                      IF (ig.GT.klon)  THEN
+                          WRITE(lunout,*) 'Attention dans readaerosol', &
+                             name_aero(id_aero), 'ig > klon', ig, klon
+                      ENDIF
+
+                      var(ig,k,it,id_aero) =  var_1(i,jjm+1+1-j,klev+1-k,it)
+
+                    ENDDO
+                  ENDDO
+                  IF (ig.NE.klon-1) THEN
+                      PRINT *,'Error in readaerosol', name_aero(id_aero), 'conversion'
+                      CALL abort_gcm('readaerosol','Error in readaerosol 1',1)
+                  ENDIF
+                ENDDO         ! Loop over k (vertical)
+              ENDDO           ! Loop over it (months)
+              
+          ENDIF               ! Had to read new data?
+           
+      
+          ! Interpolate to actual day:
+          IF (iday.LT.im*30-15) THEN         
+              ! in the first half of the month use month before and actual month
+              im2=im-1
+              day2 = im2*30-15
+              day1 = im2*30+15
+              IF (im2.LE.0) THEN 
+                  ! the month is january, thus the month before december
+                  im2=12
+              ENDIF
+              DO k=1,klev
+                DO i=1,klon
+                  massvar(i,k) = &
+                     var(i,k,im2,id_aero)- FLOAT(iday-day2)/FLOAT(day1-day2) * &
+                     (var(i,k,im2,id_aero) - var(i,k,im,id_aero))
+
+                  IF (massvar(i,k).LT.0.) THEN
+                      IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2
+                      IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN
+                          PRINT*, name_aero(id_aero),'(i,k,im2)- ', &
+                             name_aero(id_aero),'(i,k,im)',         &
+                             var(i,k,im2,id_aero) - var(i,k,im,id_aero)
+                      ENDIF
+
+                      IF (day1-day2.LT.0.) WRITE(lunout,*) 'day1-day2',day1-day2
+                      WRITE(lunout,*) 'stop',name_aero(id_aero)
+                      CALL abort_gcm('readaerosol','Error in interpolation 2',1)
+                  ENDIF
+                ENDDO
+              ENDDO
+          ELSE 
+              ! the second half of the month
+              im2=im+1
+              IF (im2.GT.12) THEN
+                  ! the month is december, the following thus january
+                  im2=1
+              ENDIF
+              day2 = im*30-15
+              day1 = im*30+15
+              DO k=1,klev
+                DO i=1,klon 
+                  massvar(i,k) = &
+                     var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
+                     (var(i,k,im2,id_aero) - var(i,k,im,id_aero))
+
+                  IF (massvar(i,k).LT.0.) THEN
+                      IF (iday-day2.LT.0.) &
+                         WRITE(lunout,*) 'iday-day2',iday-day2 
+                      IF (var(i,k,im2,id_aero) -  var(i,k,im,id_aero).LT.0.) THEN
+                          WRITE(lunout,*) name_aero(id_aero),'(i,k,im2)-', &
+                             name_aero(id_aero),'(i,k,im)',           &
+                             var(i,k,im2,id_aero) - var(i,k,im,id_aero)
+                      ENDIF
+                      IF (day1-day2.LT.0.) &
+                         WRITE(lunout,*) 'day1-day2',day1-day2
+                      WRITE(lunout,*) 'stop',name_aero(id_aero)
+                      CALL abort_gcm('readaerosol','Error in interpolation 3',1)
+                  ENDIF
+                ENDDO
+              ENDDO
+          ENDIF
+          
+      
+          ! The massvar concentration [molec cm-3] is read in. 
+          ! Convert it into mass [ug VAR/m3]
+          ! masse_var in [g/mol], n_avogadro in [molec/mol]
+          ! The sulfate mass [ug VAR/m3] is read in. 
+          DO k=1,klev 
+            DO i=1,klon 
+              var_out(i,k,id_aero) = massvar(i,k) 
+              IF (var_out(i,k,id_aero).LT.0) THEN 
+                  PRINT*, 'Attention il y a un probleme dans readaerosol', &
+                     name_aero(id_aero), 'la masse est negative'
+                  CALL abort_gcm('readaerosol','Error in readaerosol 4',1)
+              ENDIF
+            ENDDO
+          ENDDO
+      ELSE                    ! if no new day, use old data:
+          DO k=1,klev 
+            DO i=1,klon 
+              massvar(i,k) = var_out(i,k,id_aero) 
+              IF (var_out(i,k,id_aero).LT.0)   THEN 
+                  PRINT*, 'Attention il y a un probleme dans readaerosol', &
+                     name_aero(id_aero), 'la masse est negative'
+                  CALL abort_gcm('readaerosol','Error in readaerosol 5',1) 
+              ENDIF
+            ENDDO
+          ENDDO
+          
+           
+      ENDIF                   ! Did I have to do anything (was it a new day?)
+      
+  ENDIF                      ! phy_rank==0
+  
+!$OMP END MASTER 
+  CALL Scatter(massvar,massvar_p)             
+
+END SUBROUTINE readaerosol
+      
+      
+      
+      
+      
+!-----------------------------------------------------------------------------
+! Read in /calculate pre-industrial values of sulfate      
+!-----------------------------------------------------------------------------
+      
+SUBROUTINE readaerosol_preind (id_aero, r_day, first, pi_massvar_p)
+  
+  USE dimphy, ONLY : klev
+  USE mod_grid_phy_lmdz, klon=>klon_glo
+  USE mod_phys_lmdz_para
+  IMPLICIT NONE
+      
+  ! Content: 
+  ! --------
+  ! This routine reads in monthly mean values of massvar aerosols and 
+  ! returns a linearly interpolated dayly-mean field.      
+  ! 
+  ! It does so for the preindustriel values of the massvar, to a large part
+  ! analogous to the routine readmassvar above.      
+  !
+  ! Only Pb: Variables must be saved and don t have to be overwritten!
+  !      
+  ! Author:
+  ! -------
+  ! Celine Deandreis & Anne Cozic (LSCE) 2009
+  ! Johannes Quaas (quaas@lmd.jussieu.fr) 
+  ! 26/06/01
+  !
+  ! Modifications:
+  ! --------------
+  ! see above 
+
+      
+  ! Include-files:
+  ! --------------     
+
+  INCLUDE "YOMCST.h"
+  INCLUDE "chem.h"      
+  INCLUDE "dimensions.h"      
+  INCLUDE "temps.h"      
+  INCLUDE "iniprint.h"
+ 
+  ! Input:
+  ! ------
+  INTEGER, INTENT(in) ::  id_aero
+  REAL, INTENT(in)    ::  r_day          ! Day of integration
+  LOGICAL, INTENT(in) ::  first          ! First timestep (and therefore initialization necessary)
+
+
+  !      
+  ! Output:      
+  ! -------     
+  REAL, DIMENSION(klon_omp,klev), INTENT(out) ::  pi_massvar_p ! Number conc. massvar (monthly mean data, from file)
+  
+
+  !     
+  ! Local Variables:
+  ! ----------------      
+  INTEGER                                      :: i, ig, k, it
+  INTEGER                                      :: j, iday, iyr
+  INTEGER, PARAMETER                           :: ny=jjm+1
+  INTEGER                                      :: im, day1, day2, im2
+
+  REAL, DIMENSION(iim, jjm+1, klev, 12)        :: pi_var_1
+  REAL, DIMENSION(klon,klev)                   :: pi_massvar
+  REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE  :: pi_var ! VAR in right dimension
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE    :: pi_var_out
+!$OMP THREADPRIVATE(pi_var,pi_var_out)            
+
+  CHARACTER(len=4)                :: cyear
+  CHARACTER(len=7),DIMENSION(8)   :: name_aero
+  LOGICAL                         :: lnewday
+  LOGICAL,SAVE                    :: first2=.TRUE.
+!$OMP THREADPRIVATE(first2)
+  
+
+! Variable pour definir a partir d'un numero d'aerosol, son nom
+  name_aero(1) = "SSSSM  "
+  name_aero(2) = "ASSSM  "
+  name_aero(3) = "ASBCM  "
+  name_aero(4) = "ASPOMM "
+  name_aero(5) = "SO4    "
+  name_aero(6) = "CIDUSTM"
+  name_aero(7) = "AIBCM  "
+  name_aero(8) = "AIPOMM "
+
+
+!$OMP MASTER
+  IF (first2) THEN
+      ALLOCATE( pi_var(klon, klev, 12,8) )
+      ALLOCATE( pi_var_out(klon, klev,8))
+      first2=.FALSE.
+  ENDIF
+
+  IF (is_mpi_root) THEN
+      iday = INT(r_day) 
+      ! Get the year of the run
+      iyr  = iday/360
+      ! Get the day of the actual year:
+      iday = iday-iyr*360
+      ! 0.02 is about 0.5/24, namly less than half an hour
+      lnewday = (r_day-FLOAT(iday).LT.0.02)
+      
+      !     ---------------------------------------------
+      !     All has to be done only, if a new day begins!       
+      !     ---------------------------------------------
+      
+      IF (lnewday.OR.first) THEN
+          
+          im = iday/30 +1     ! the actual month
+          ! annee_ref is the initial year (defined in temps.h)
+          iyr = iyr + annee_ref      
+
+          IF (first) THEN
+              cyear='.nat'
+              CALL get_aero_fromfile(cyear,pi_var_1, name_aero(id_aero))
+              
+              ! Transform the horizontal 2D-field into the physics-field
+              ! (Also the levels and the latitudes have to be inversed)
+              
+              ! Initialize field
+              DO it=1,12
+                DO k=1,klev
+                  DO i=1,klon
+                    pi_var(i,k,it,id_aero)=0.
+                  ENDDO
+                ENDDO
+              ENDDO
+               
+              WRITE(lunout,*) 'preind: finished reading', FLOAT(iim)
+              DO it=1,12
+                DO k=1, klev         
+                  ! a) at the poles, use the zonal mean:
+                  DO i=1,iim
+                    ! North pole
+                    pi_var(1,k,it,id_aero)    = &
+                       pi_var(1,k,it,id_aero) + pi_var_1(i,jjm+1,klev+1-k,it)
+                    ! South pole
+                    pi_var(klon,k,it,id_aero) = &
+                       pi_var(klon,k,it,id_aero) + pi_var_1(i,1,klev+1-k,it)
+                  ENDDO
+                  pi_var(1,k,it,id_aero)    = pi_var(1,k,it,id_aero)/FLOAT(iim)
+                  pi_var(klon,k,it,id_aero) = pi_var(klon,k,it,id_aero)/FLOAT(iim)
+                  
+                  ! b) the values between the poles:
+                  ig=1
+                  DO j=2,jjm
+                    DO i=1,iim
+                      ig=ig+1
+                      IF (ig.GT.klon)  THEN
+                          WRITE(lunout,*) 'Attention dans readaerosol_preind', &
+                             name_aero(id_aero), 'ig > klon', ig, klon
+                      ENDIF
+                      pi_var(ig,k,it,id_aero) = & 
+                         pi_var_1(i,jjm+1+1-j,klev+1-k,it)
+                    ENDDO
+                  ENDDO
+                  IF (ig.NE.klon-1) THEN
+                      WRITE(lunout,*) 'Error in readaerosol_preind (', name_aero(id_aero), ' conversion)'
+                      CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 1',1) 
+                  ENDIF
+                ENDDO         ! Loop over k (vertical)
+              ENDDO             ! Loop over it (months)
+               
+          ENDIF                ! Had to read new data?
+            
+            
+                                ! Interpolate to actual day:
+          IF (iday.LT.im*30-15) THEN          
+              ! in the first half of the month use month before and actual month
+              im2=im-1 
+              day1 = im2*30+15 
+              day2 = im2*30-15 
+              IF (im2.LE.0) THEN  
+                  ! the month is january, thus the month before december
+                  im2=12 
+              ENDIF
+              DO k=1,klev 
+                DO i=1,klon
+                  pi_massvar(i,k) = &
+                     pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
+                     (pi_var(i,k,im2,id_aero) -  pi_var(i,k,im,id_aero)) 
+
+                  IF (pi_massvar(i,k).LT.0.) THEN 
+                      IF (iday-day2.LT.0.)  WRITE(lunout,*) 'iday-day2',iday-day2 
+                      IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN
+                          WRITE(lunout,*) 'pi_',name_aero(id_aero),'(i,k,im2) - pi_', &
+                             name_aero(id_aero),'(i,k,im)', &
+                             pi_var(i,k,im2,id_aero) -pi_var(i,k,im,id_aero) 
+                      ENDIF
+                      IF (day1-day2.LT.0.)WRITE(lunout,*) 'day1-day2',day1-day2 
+                      PRINT *, 'pi_',name_aero(id_aero)
+                      CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 2',1)
+                  ENDIF
+                ENDDO
+              ENDDO
+          ELSE
+
+              ! the second half of the month
+              im2=im+1
+              day1 = im*30+15
+              IF (im2.GT.12) THEN
+                  ! the month is december, the following thus january
+                  im2=1
+              ENDIF
+              day2 = im*30-15
+              
+              DO k=1,klev
+                DO i=1,klon
+                  pi_massvar(i,k) = &
+                     pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
+                     (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero))
+                  IF (pi_massvar(i,k).LT.0.) THEN
+                      IF (iday-day2.LT.0.)  WRITE(lunout,*) 'iday-day2',iday-day2 
+                      IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN
+                          WRITE(lunout,*) 'pi_', name_aero(id_aero),'(i,k,im2) - pi_',&
+                             name_aero(id_aero), '(i,k,im)',&
+                             pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)
+                      ENDIF
+                      IF (day1-day2.LT.0.) &
+                         WRITE(lunout,*) 'day1-day2',day1-day2
+                      PRINT *, 'pi_',name_aero(id_aero)
+                      CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 3',1)
+                  ENDIF
+                ENDDO
+              ENDDO
+          ENDIF
+          
+         
+          ! The massvar concentration [molec cm-3] is read in. 
+          ! Convert it into mass [ug VAR/m3]
+          ! masse_var in [g/mol], n_avogadro in [molec/mol]
+          DO k=1,klev
+            DO i=1,klon
+              pi_var_out(i,k,id_aero) = pi_massvar(i,k)
+            ENDDO
+          ENDDO
+          
+      ELSE                   ! If no new day, use old data:
+          DO k=1,klev
+            DO i=1,klon
+              pi_massvar(i,k) = pi_var_out(i,k,id_aero)            
+            ENDDO
+          ENDDO
+      ENDIF                  ! Was this the beginning of a new day?
+   ENDIF                     ! is_mpi_root
+      
+!$OMP END MASTER
+   CALL Scatter(pi_massvar,pi_massvar_p)            
+      
+ END SUBROUTINE readaerosol_preind
+      
+      
+      
+!-----------------------------------------------------------------------------
+! Routine for reading VAR data from files
+!-----------------------------------------------------------------------------
+            
+
+SUBROUTINE get_aero_fromfile (cyr, var, varname)
+  USE dimphy
+  IMPLICIT NONE
+      
+  INCLUDE "netcdf.inc"
+  INCLUDE "dimensions.h"      
+  INCLUDE "iniprint.h"
+
+  CHARACTER(len=4), INTENT(in)                       ::  cyr
+  REAL, DIMENSION(iim, jjm+1, klev, 12), INTENT(out) ::  var
+  CHARACTER(len=*), INTENT(in)                       ::  varname
+  
+
+  CHARACTER(len=30)     :: fname
+  CHARACTER(len=30)     :: cvar
+  INTEGER, DIMENSION(3) :: START, COUNT
+  INTEGER               :: STATUS, NCID, VARID
+  INTEGER               :: imth, i, j, k
+  INTEGER, PARAMETER    :: ny=jjm+1
+  REAL, DIMENSION(iim, ny, klev) :: varmth
+!
+!
+!  
+  fname = TRIM(varname)//'.run'//cyr//'.cdf'
+  
+  WRITE(lunout,*) 'reading ', TRIM(fname)
+  STATUS = NF_OPEN (TRIM(fname), NF_NOWRITE, NCID)
+  IF (STATUS .NE. NF_NOERR) THEN
+      WRITE(lunout,*) 'err in open ',status
+      CALL abort_gcm('get_aero_fromfile','Error in opening file',1)
+  ENDIF
+  DO imth=1, 12
+    IF (imth.EQ.1) THEN
+        cvar=TRIM(varname)//'JAN'
+    ELSEIF (imth.EQ.2) THEN
+        cvar=TRIM(varname)//'FEB'
+    ELSEIF (imth.EQ.3) THEN
+        cvar=TRIM(varname)//'MAR'
+    ELSEIF (imth.EQ.4) THEN
+        cvar=TRIM(varname)//'APR'
+    ELSEIF (imth.EQ.5) THEN
+        cvar=TRIM(varname)//'MAY'
+    ELSEIF (imth.EQ.6) THEN
+        cvar=TRIM(varname)//'JUN'
+    ELSEIF (imth.EQ.7) THEN
+        cvar=TRIM(varname)//'JUL'
+    ELSEIF (imth.EQ.8) THEN
+        cvar=TRIM(varname)//'AUG'
+    ELSEIF (imth.EQ.9) THEN
+        cvar=TRIM(varname)//'SEP'
+    ELSEIF (imth.EQ.10) THEN
+        cvar=TRIM(varname)//'OCT'
+    ELSEIF (imth.EQ.11) THEN
+        cvar=TRIM(varname)//'NOV'
+    ELSEIF (imth.EQ.12) THEN
+        cvar=TRIM(varname)//'DEC'
+    ENDIF
+    start(1)=1
+    start(2)=1
+    start(3)=1
+    COUNT(1)=iim
+    COUNT(2)=ny
+    COUNT(3)=klev
+    STATUS = NF_INQ_VARID (NCID, TRIM(cvar), VARID)
+    WRITE(lunout,*) ncid,imth,TRIM(cvar), varid
+    
+    IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read ',status      
+
+#ifdef NC_DOUBLE
+    status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT,varmth)
+#else
+    status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, varmth)
+#endif
+    IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read data',status
+         
+    DO k=1,klev
+      DO j=1,jjm+1
+        DO i=1,iim
+          IF (varmth(i,j,k).LT.0.) THEN
+              WRITE(lunout,*) 'this is shit'
+              WRITE(lunout,*) varname,'(',i,j,k,') =',varmth(i,j,k)
+          ENDIF
+          var(i,j,k,imth)=varmth(i,j,k)
+        ENDDO
+      ENDDO
+    ENDDO
+  ENDDO
+  
+  STATUS = NF_CLOSE(NCID)
+  IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in closing file',status
+  
+
+END SUBROUTINE get_aero_fromfile
+      
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readsulfate.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readsulfate.F	(revision 1145)
+++ 	(revision )
@@ -1,660 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE readsulfate (r_day, first, sulfate_p)
-      USE dimphy, ONLY : klev
-      USE mod_grid_phy_lmdz, klon=>klon_glo
-      USE mod_phys_lmdz_para
-      IMPLICIT none
-      
-c Content: 
-c --------
-c This routine reads in monthly mean values of sulfate aerosols and 
-c returns a linearly interpolated dayly-mean field.      
-c 
-c
-c Author:
-c -------
-c Johannes Quaas (quaas@lmd.jussieu.fr) 
-c 26/04/01
-c
-c Modifications:
-c --------------
-c 21/06/01: Make integrations of more than one year possible ;-)     
-c           ATTENTION!! runs are supposed to start with Jan, 1. 1930
-c                       (rday=1)      
-c
-c 27/06/01: Correction: The model always has 360 days per year!
-c 27/06/01: SO4 concentration rather than mixing ratio      
-c 27/06/01: 10yr-mean-values to interpolate     
-c 20/08/01: Correct the error through integer-values in interpolations      
-c 21/08/01: Introduce flag to read in just one decade
-c      
-c Include-files:
-c --------------     
-#include "YOMCST.h"
-#include "chem.h"      
-#include "dimensions.h"      
-#include "temps.h"      
-#include "clesphys.h"
-#include "iniprint.h"
-c 
-c Input:
-c ------
-      REAL  r_day                   ! Day of integration
-      LOGICAL first                 ! First timestep 
-                                    ! (and therefore initialization necessary)
-c      
-c Output:      
-c -------     
-      REAL  sulfate_p(klon_omp,klev)
-      REAL  sulfate (klon, klev)  ! Mass of sulfate (monthly mean data, 
-                                  !  from file) [ug SO4/m3]
-c      
-c Local Variables:
-c ----------------      
-      INTEGER i, ig, k, it
-      INTEGER j, iday, ny, iyr, iyr1, iyr2
-      parameter (ny=jjm+1)
-      
-CJLD      INTEGER idec1, idec2 ! The two decadal data read ini
-      CHARACTER*4 cyear
-      
-      INTEGER im, day1, day2, im2
-      REAL so4_1(iim, jjm+1, klev, 12)
-      REAL so4_2(iim, jjm+1, klev, 12)   ! The sulfate distributions
-      
-      REAL, allocatable,save :: so4(:, :, :)  ! SO4 in right dimension
-      REAL, allocatable,save :: so4_out(:, :)
-c$OMP THREADPRIVATE(so4,so4_out)
-      
-      LOGICAL lnewday 
-      LOGICAL lonlyone
-      PARAMETER (lonlyone=.FALSE.)
-      logical,save :: first2=.true.
-c$OMP THREADPRIVATE(first2)
-
-c$OMP MASTER
-      if (first2) then
-     
-        allocate( so4(klon, klev, 12) )
-        allocate( so4_out(klon, klev))
-        first2=.false.
-	
-      endif
-
-      if (is_mpi_root) then
-
-        IF (aer_type /= 'actuel  ' .AND. aer_type /= 'preind  ' .AND.   &
-     &      aer_type /= 'scenario') THEN
-          WRITE(lunout,*)' *** Warning ***'
-          WRITE(lunout,*)'Option aer_type pour les aerosols = ',        &
-     &        aer_type
-          WRITE(lunout,*)'Cas non prevu, force a preind'
-          aer_type = 'preind  '
-        ENDIF
-            
-      iday = INT(r_day) 
-      
-      ! Get the year of the run
-      iyr  = iday/360
-      
-      ! Get the day of the actual year:
-      iday = iday-iyr*360
-      
-      ! 0.02 is about 0.5/24, namly less than half an hour
-      lnewday = (r_day-FLOAT(iday).LT.0.02)
-      
-! ---------------------------------------------
-! All has to be done only, if a new day begins!       
-! ---------------------------------------------
-
-      IF (lnewday.OR.first) THEN
-         
-      im = iday/30 +1 ! the actual month
-      ! annee_ref is the initial year (defined in temps.h)
-      iyr = iyr + annee_ref
-      
-      ! Do I have to read new data? (Is this the first day of a year?)
-      IF (first.OR.iday.EQ.1.) THEN 
-      ! Initialize values
-      DO it=1,12
-      DO k=1,klev
-         DO i=1,klon
-            so4(i,k,it)=0.
-         ENDDO
-      ENDDO
-      ENDDO
-
-
-
-      IF (aer_type == 'actuel  ') then
-        cyear='1980'
-        CALL getso4fromfile(cyear, so4_1)
-      ELSE IF (aer_type == 'preind  ') THEN
-        cyear='.nat'
-        CALL getso4fromfile(cyear, so4_1)
-      ELSE
-        IF (iyr .lt. 1850) THEN
-           cyear='.nat'
-           WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
-           CALL getso4fromfile(cyear, so4_1)
-        ELSE IF (iyr .ge. 2100) THEN
-           cyear='2100'
-           WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
-           CALL getso4fromfile(cyear, so4_1)
-        ELSE
-
-      	! Read in data:
-        ! a) from actual 10-yr-period
-
-          IF (iyr.LT.1900) THEN
-             iyr1 = 1850
-             iyr2 = 1900
-          ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
-             iyr1 = 1900
-             iyr2 = 1920
-          ELSE 
-             iyr1 = INT(iyr/10)*10
-             iyr2 = INT(1+iyr/10)*10
-          ENDIF
-          WRITE(cyear,'(I4)') iyr1
-        ENDIF
-        WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
-        CALL getso4fromfile(cyear, so4_1)
-
-      
-      ! If to read two decades:
-        IF (.NOT.lonlyone) THEN
-         
-      ! b) from the next following one
-          WRITE(cyear,'(I4)') iyr2
-          WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
-          CALL getso4fromfile(cyear, so4_2)
-
- 
-      ! Interpolate linarily to the actual year:
-        DO it=1,12
-           DO k=1,klev
-              DO j=1,jjm
-                 DO i=1,iim
-                    so4_1(i,j,k,it)=so4_1(i,j,k,it)
-     .                 - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1)
-     .                 * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
-                 ENDDO
-              ENDDO
-           ENDDO
-        ENDDO                           
-
-
-        ENDIF !lonlyone    
-      ENDIF !aer_type
-      
-      ! Transform the horizontal 2D-field into the physics-field
-      ! (Also the levels and the latitudes have to be inversed)
-      
-      DO it=1,12
-      DO k=1, klev         
-         ! a) at the poles, use the zonal mean:
-         DO i=1,iim
-            ! North pole
-            so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
-            ! South pole
-            so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
-         ENDDO
-         so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
-         so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)
-      
-         ! b) the values between the poles:
-         ig=1
-         DO j=2,jjm
-            DO i=1,iim
-               ig=ig+1
-               if (ig.gt.klon) write (*,*) 'shit'
-               so4(ig,k,it) = so4_1(i,jjm+1+1-j,klev+1-k,it)
-            ENDDO
-         ENDDO
-         IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
-      ENDDO ! Loop over k (vertical)
-      ENDDO ! Loop over it (months)
-               
-
-      ENDIF ! Had to read new data?
-      
-      
-      ! Interpolate to actual day:
-      IF (iday.LT.im*30-15) THEN         
-         ! in the first half of the month use month before and actual month
-         im2=im-1
-         day2 = im2*30-15
-         day1 = im2*30+15
-         IF (im2.LE.0) THEN 
-            ! the month is january, thus the month before december
-            im2=12
-         ENDIF
-         DO k=1,klev
-            DO i=1,klon
-               sulfate(i,k) = so4(i,k,im2)  
-     .              - FLOAT(iday-day2)/FLOAT(day1-day2)
-     .              * (so4(i,k,im2) - so4(i,k,im))
-               IF (sulfate(i,k).LT.0.) THEN
-                  IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
-                  IF (so4(i,k,im2) - so4(i,k,im).LT.0.)
-     . write(*,*) 'so4(i,k,im2) - so4(i,k,im)',
-     . so4(i,k,im2) - so4(i,k,im)
-                  IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
-                  stop 'sulfate'
-               endif
-            ENDDO
-         ENDDO
-      ELSE 
-         ! the second half of the month
-         im2=im+1
-         IF (im2.GT.12) THEN
-            ! the month is december, the following thus january
-            im2=1
-         ENDIF
-         day2 = im*30-15
-         day1 = im*30+15
-         DO k=1,klev
-            DO i=1,klon
-               sulfate(i,k) = so4(i,k,im2)  
-     .              - FLOAT(iday-day2)/FLOAT(day1-day2)
-     .              * (so4(i,k,im2) - so4(i,k,im))
-               IF (sulfate(i,k).LT.0.) THEN
-                  IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
-                  IF (so4(i,k,im2) - so4(i,k,im).LT.0.)
-     . write(*,*) 'so4(i,k,im2) - so4(i,k,im)',
-     . so4(i,k,im2) - so4(i,k,im)
-                  IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
-                  stop 'sulfate'
-               endif
-            ENDDO
-         ENDDO
-      ENDIF
-
-      
-CJLD      ! The sulfate concentration [molec cm-3] is read in. 
-CJLD      ! Convert it into mass [ug SO4/m3]
-CJLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
-      ! The sulfate mass [ug SO4/m3] is read in. 
-      DO k=1,klev
-         DO i=1,klon
-CJLD            sulfate(i,k) = sulfate(i,k)*masse_so4
-CJLD     .           /n_avogadro*1.e+12
-            so4_out(i,k) = sulfate(i,k)
-            IF (so4_out(i,k).LT.0) 
-     .          stop 'WAS SOLL DER SCHEISS ? '
-         ENDDO
-      ENDDO
-      ELSE ! if no new day, use old data:
-      DO k=1,klev
-         DO i=1,klon
-            sulfate(i,k) = so4_out(i,k)
-            IF (so4_out(i,k).LT.0) 
-     .          stop 'WAS SOLL DER SCHEISS ? '
-         ENDDO
-      ENDDO
-         
-
-      ENDIF ! Did I have to do anything (was it a new day?)
-
-      endif   ! phy_rank==0
-      
-c$OMP END MASTER
-      call Scatter(sulfate,sulfate_p)            
-
-      RETURN
-      END
-
-      
-      
-      
-      
-c-----------------------------------------------------------------------------
-c Read in /calculate pre-industrial values of sulfate      
-c-----------------------------------------------------------------------------
-      
-      SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate_p)
-      USE dimphy, ONLY : klev
-      USE mod_grid_phy_lmdz, klon=>klon_glo
-      USE mod_phys_lmdz_para
-      IMPLICIT none
-      
-c Content: 
-c --------
-c This routine reads in monthly mean values of sulfate aerosols and 
-c returns a linearly interpolated dayly-mean field.      
-c 
-c It does so for the preindustriel values of the sulfate, to a large part
-c analogous to the routine readsulfate above.      
-c
-c Only Pb: Variables must be saved and don t have to be overwritten!
-c      
-c Author:
-c -------
-c Johannes Quaas (quaas@lmd.jussieu.fr) 
-c 26/06/01
-c
-c Modifications:
-c --------------
-c see above 
-c      
-c Include-files:
-c --------------     
-#include "YOMCST.h"
-#include "chem.h"      
-#include "dimensions.h"      
-#include "temps.h"      
-c 
-c Input:
-c ------
-      REAL  r_day                   ! Day of integration
-      LOGICAL first                 ! First timestep 
-                                    ! (and therefore initialization necessary)
-c      
-c Output:      
-c -------     
-      REAL  pi_sulfate_p (klon_omp, klev)  
-                                 
-      REAL  pi_sulfate (klon, klev)  ! Number conc. sulfate (monthly mean data, 
-                                  !  from fil
-c      
-c Local Variables:
-c ----------------      
-      INTEGER i, ig, k, it
-      INTEGER j, iday, ny, iyr
-      parameter (ny=jjm+1)
-      
-      INTEGER im, day1, day2, im2
-      REAL pi_so4_1(iim, jjm+1, klev, 12)
-
-      REAL, allocatable,save :: pi_so4(:, :, :)  ! SO4 in right dimension
-      REAL, allocatable,save :: pi_so4_out(:, :)
-c$OMP THREADPRIVATE(pi_so4,pi_so4_out)            
-      
-      CHARACTER*4 cyear
-      LOGICAL lnewday
-      logical,save :: first2=.true.
-c$OMP THREADPRIVATE(first2)
-
-c$OMP MASTER
-      if (first2) then
-     
-        allocate( pi_so4(klon, klev, 12) )
-        allocate( pi_so4_out(klon, klev))
-        first2=.false.
-	
-      endif
-
-      if (is_mpi_root) then
-    
-      
-
-      iday = INT(r_day) 
-      
-      ! Get the year of the run
-      iyr  = iday/360
-      
-      ! Get the day of the actual year:
-      iday = iday-iyr*360
-      
-      ! 0.02 is about 0.5/24, namly less than half an hour
-      lnewday = (r_day-FLOAT(iday).LT.0.02)
-      
-! ---------------------------------------------
-! All has to be done only, if a new day begins!       
-! ---------------------------------------------
-
-      IF (lnewday.OR.first) THEN
-         
-      im = iday/30 +1 ! the actual month
-      
-      ! annee_ref is the initial year (defined in temps.h)
-      iyr = iyr + annee_ref      
-      
-      
-      IF (first) THEN
-         cyear='.nat'
-         CALL getso4fromfile(cyear,pi_so4_1)
-
-               ! Transform the horizontal 2D-field into the physics-field
-               ! (Also the levels and the latitudes have to be inversed)
-
-         ! Initialize field
-         DO it=1,12
-            DO k=1,klev
-               DO i=1,klon
-                  pi_so4(i,k,it)=0.
-               ENDDO
-            ENDDO
-         ENDDO
-         
-         write (*,*) 'preind: finished reading', FLOAT(iim)
-      DO it=1,12
-      DO k=1, klev         
-         ! a) at the poles, use the zonal mean:
-         DO i=1,iim
-            ! North pole
-            pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it)
-            ! South pole
-           pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it)
-         ENDDO
-         pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim)
-         pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim)
-      
-         ! b) the values between the poles:
-         ig=1
-         DO j=2,jjm
-            DO i=1,iim
-               ig=ig+1
-               if (ig.gt.klon) write (*,*) 'shit'
-               pi_so4(ig,k,it) = pi_so4_1(i,jjm+1+1-j,klev+1-k,it)
-            ENDDO
-         ENDDO
-         IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
-      ENDDO ! Loop over k (vertical)
-      ENDDO ! Loop over it (months)
-
-      ENDIF                     ! Had to read new data?
-      
-      
-      ! Interpolate to actual day:
-      IF (iday.LT.im*30-15) THEN         
-         ! in the first half of the month use month before and actual month
-         im2=im-1
-         day1 = im2*30+15
-         day2 = im2*30-15
-         IF (im2.LE.0) THEN 
-            ! the month is january, thus the month before december
-            im2=12
-         ENDIF
-         DO k=1,klev
-            DO i=1,klon
-               pi_sulfate(i,k) = pi_so4(i,k,im2)  
-     .              - FLOAT(iday-day2)/FLOAT(day1-day2)
-     .              * (pi_so4(i,k,im2) - pi_so4(i,k,im))
-               IF (pi_sulfate(i,k).LT.0.) THEN
-                  IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
-                  IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.)
-     . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)',
-     . pi_so4(i,k,im2) - pi_so4(i,k,im)
-                  IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
-                  stop 'pi_sulfate'
-               endif
-            ENDDO
-         ENDDO
-      ELSE 
-         ! the second half of the month
-         im2=im+1
-         day1 = im*30+15
-         IF (im2.GT.12) THEN
-            ! the month is december, the following thus january
-            im2=1
-         ENDIF
-         day2 = im*30-15
-         
-         DO k=1,klev
-            DO i=1,klon
-               pi_sulfate(i,k) = pi_so4(i,k,im2)  
-     .              - FLOAT(iday-day2)/FLOAT(day1-day2)
-     .              * (pi_so4(i,k,im2) - pi_so4(i,k,im))
-               IF (pi_sulfate(i,k).LT.0.) THEN
-                  IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
-                  IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.)
-     . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)',
-     . pi_so4(i,k,im2) - pi_so4(i,k,im)
-                  IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
-                  stop 'pi_sulfate'
-               endif
-            ENDDO
-         ENDDO
-      ENDIF
-
-      
-CJLD      ! The sulfate concentration [molec cm-3] is read in. 
-CJLD      ! Convert it into mass [ug SO4/m3]
-CJLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
-      DO k=1,klev
-         DO i=1,klon
-CJLD            pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
-CJLD     .           /n_avogadro*1.e+12
-            pi_so4_out(i,k) = pi_sulfate(i,k)
-         ENDDO
-      ENDDO
-      
-      ELSE ! If no new day, use old data:
-      DO k=1,klev
-         DO i=1,klon
-            pi_sulfate(i,k) = pi_so4_out(i,k)            
-         ENDDO
-      ENDDO
-         
-
-      ENDIF ! Was this the beginning of a new day?
-
-      endif   ! is_mpi_root==0
-      
-c$OMP END MASTER
-      call Scatter(pi_sulfate,pi_sulfate_p)            
-
-      RETURN
-      END
-
-      
-      
-      
-      
-      
-      
-      
-      
-      
-c-----------------------------------------------------------------------------
-c Routine for reading SO4 data from files
-c-----------------------------------------------------------------------------
-            
-
-      SUBROUTINE getso4fromfile (cyr, so4)
-      USE dimphy
-#include "netcdf.inc"
-#include "dimensions.h"      
-      CHARACTER*15 fname
-      CHARACTER*4 cyr
-      
-      CHARACTER*6 cvar
-      INTEGER START(3), COUNT(3)
-      INTEGER  STATUS, NCID, VARID
-      INTEGER imth, i, j, k, ny
-      PARAMETER (ny=jjm+1)
-      
-            
-      REAL so4mth(iim, ny, klev)
-      REAL so4(iim, ny, klev, 12)
-
- 
-      fname = 'so4.run'//cyr//'.cdf'
-
-      write (*,*) 'reading ', fname
-      STATUS = NF_OPEN (fname, NF_NOWRITE, NCID)
-      IF (STATUS .NE. NF_NOERR) write (*,*) 'err in open ',status
-            
-      DO imth=1, 12
-         IF (imth.eq.1) THEN
-            cvar='SO4JAN'
-         ELSEIF (imth.eq.2) THEN
-            cvar='SO4FEB'
-         ELSEIF (imth.eq.3) THEN
-            cvar='SO4MAR'
-         ELSEIF (imth.eq.4) THEN
-            cvar='SO4APR'
-         ELSEIF (imth.eq.5) THEN
-            cvar='SO4MAY'
-         ELSEIF (imth.eq.6) THEN
-            cvar='SO4JUN'
-         ELSEIF (imth.eq.7) THEN
-            cvar='SO4JUL'
-         ELSEIF (imth.eq.8) THEN
-            cvar='SO4AUG'
-         ELSEIF (imth.eq.9) THEN
-            cvar='SO4SEP'
-         ELSEIF (imth.eq.10) THEN
-            cvar='SO4OCT'
-         ELSEIF (imth.eq.11) THEN
-            cvar='SO4NOV'
-         ELSEIF (imth.eq.12) THEN
-            cvar='SO4DEC'
-         ENDIF
-         start(1)=1
-         start(2)=1
-         start(3)=1
-         count(1)=iim
-         count(2)=ny
-         count(3)=klev
-c         write(*,*) 'here i am'
-         STATUS = NF_INQ_VARID (NCID, cvar, VARID)
-         write (*,*) ncid,imth,cvar, varid
-
-         IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read ',status      
-
-#ifdef NC_DOUBLE
-         status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT, so4mth)
-#else
-         status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, so4mth)
-#endif
-         IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read data',status
-         
-         DO k=1,klev
-            DO j=1,jjm+1
-               DO i=1,iim
-                  IF (so4mth(i,j,k).LT.0.) then
-                     write(*,*) 'this is shit'
-                     write(*,*) 'so4(',i,j,k,') =',so4mth(i,j,k)
-                  endif
-                  so4(i,j,k,imth)=so4mth(i,j,k)
-               ENDDO
-            ENDDO
-         ENDDO
-      ENDDO
-      
-      STATUS = NF_CLOSE(NCID)
-      IF (STATUS .NE. NF_NOERR) write (*,*) 'err in closing file',status
-
-
-      END ! subroutine getso4fromfile
-      
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aero.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aero.F90	(revision 1150)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aero.F90	(revision 1150)
@@ -0,0 +1,750 @@
+SUBROUTINE SW_AERO(PSCT, PRMU0, PFRAC, &
+     PPMB, PDP, &
+     PPSOL, PALBD, PALBP,&
+     PTAVE, PWV, PQS, POZON, PAER,&
+     PCLDSW, PTAU, POMEGA, PCG,&
+     PHEAT, PHEAT0,&
+     PALBPLA,PTOPSW,PSOLSW,PTOPSW0,PSOLSW0,&
+     ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
+     tauaero, pizaero, cgaero,&
+     PTAUA, POMEGAA,&
+     PTOPSWADAERO,PSOLSWADAERO,&
+     PTOPSWAD0AERO,PSOLSWAD0AERO,&
+     PTOPSWAIAERO,PSOLSWAIAERO,&
+     PTOPSWAERO,PTOPSW0AERO,&
+     PSOLSWAERO,PSOLSW0AERO,&
+     ok_ade, ok_aie )
+
+  USE dimphy
+  IMPLICIT NONE
+
+#include "YOMCST.h"
+  !
+  !     ------------------------------------------------------------------
+  !
+  !     PURPOSE.
+  !     --------
+  !
+  !          THIS ROUTINE COMPUTES THE SHORTWAVE RADIATION FLUXES IN TWO
+  !     SPECTRAL INTERVALS FOLLOWING FOUQUART AND BONNEL (1980).
+  !
+  !     METHOD.
+  !     -------
+  !
+  !          1. COMPUTES ABSORBER AMOUNTS                 (SWU)
+  !          2. COMPUTES FLUXES IN 1ST SPECTRAL INTERVAL  (SW1S)
+  !          3. COMPUTES FLUXES IN 2ND SPECTRAL INTERVAL  (SW2S)
+  !
+  !     REFERENCE.
+  !     ----------
+  !
+  !        SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
+  !        DOCUMENTATION, AND FOUQUART AND BONNEL (1980)
+  !
+  !     AUTHOR.
+  !     -------
+  !        JEAN-JACQUES MORCRETTE  *ECMWF*
+  !
+  !     MODIFICATIONS.
+  !     --------------
+  !        ORIGINAL : 89-07-14
+  !        95-01-01   J.-J. MORCRETTE  Direct/Diffuse Albedo
+  !        03-11-27   J. QUAAS Introduce aerosol forcings (based on BOUCHER)
+  !        09-04      A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing
+  !     ------------------------------------------------------------------
+  !
+  !* ARGUMENTS:
+  !
+  REAL*8 PSCT  ! constante solaire (valeur conseillee: 1370)
+
+  REAL*8 PPSOL(KDLON)        ! SURFACE PRESSURE (PA)
+  REAL*8 PDP(KDLON,KFLEV)    ! LAYER THICKNESS (PA)
+  REAL*8 PPMB(KDLON,KFLEV+1) ! HALF-LEVEL PRESSURE (MB)
+
+  REAL*8 PRMU0(KDLON)  ! COSINE OF ZENITHAL ANGLE
+  REAL*8 PFRAC(KDLON)  ! fraction de la journee
+
+  REAL*8 PTAVE(KDLON,KFLEV)  ! LAYER TEMPERATURE (K)
+  REAL*8 PWV(KDLON,KFLEV)    ! SPECIFI! HUMIDITY (KG/KG)
+  REAL*8 PQS(KDLON,KFLEV)    ! SATURATED WATER VAPOUR (KG/KG)
+  REAL*8 POZON(KDLON,KFLEV)  ! OZONE CONCENTRATION (KG/KG)
+  REAL*8 PAER(KDLON,KFLEV,5) ! AEROSOLS' OPTICAL THICKNESS
+
+  REAL*8 PALBD(KDLON,2)  ! albedo du sol (lumiere diffuse)
+  REAL*8 PALBP(KDLON,2)  ! albedo du sol (lumiere parallele)
+
+  REAL*8 PCLDSW(KDLON,KFLEV)    ! CLOUD FRACTION
+  REAL*8 PTAU(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS
+  REAL*8 PCG(KDLON,2,KFLEV)     ! ASYMETRY FACTOR
+  REAL*8 POMEGA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
+
+  REAL*8 PHEAT(KDLON,KFLEV) ! SHORTWAVE HEATING (K/DAY)
+  REAL*8 PHEAT0(KDLON,KFLEV)! SHORTWAVE HEATING (K/DAY) clear-sky
+  REAL*8 PALBPLA(KDLON)     ! PLANETARY ALBEDO
+  REAL*8 PTOPSW(KDLON)      ! SHORTWAVE FLUX AT T.O.A.
+  REAL*8 PSOLSW(KDLON)      ! SHORTWAVE FLUX AT SURFACE
+  REAL*8 PTOPSW0(KDLON)     ! SHORTWAVE FLUX AT T.O.A. (CLEAR-SKY)
+  REAL*8 PSOLSW0(KDLON)     ! SHORTWAVE FLUX AT SURFACE (CLEAR-SKY)
+  !
+  !* LOCAL VARIABLES:
+  !
+  REAL*8 ZOZ(KDLON,KFLEV)
+  REAL*8 ZAKI(KDLON,2)     
+  REAL*8 ZCLD(KDLON,KFLEV)
+  REAL*8 ZCLEAR(KDLON) 
+  REAL*8 ZDSIG(KDLON,KFLEV)
+  REAL*8 ZFACT(KDLON)
+  REAL*8 ZFD(KDLON,KFLEV+1)
+  REAL*8 ZFDOWN(KDLON,KFLEV+1)
+  REAL*8 ZFU(KDLON,KFLEV+1)
+  REAL*8 ZFUP(KDLON,KFLEV+1)
+  REAL*8 ZRMU(KDLON)
+  REAL*8 ZSEC(KDLON)
+  REAL*8 ZUD(KDLON,5,KFLEV+1)
+  REAL*8 ZCLDSW0(KDLON,KFLEV)
+
+  REAL*8 ZFSUP(KDLON,KFLEV+1)
+  REAL*8 ZFSDN(KDLON,KFLEV+1)
+  REAL*8 ZFSUP0(KDLON,KFLEV+1)
+  REAL*8 ZFSDN0(KDLON,KFLEV+1)
+
+  INTEGER inu, jl, jk, i, k, kpl1
+
+  INTEGER swpas  ! Every swpas steps, sw is calculated
+  PARAMETER(swpas=1)
+
+  INTEGER itapsw
+  LOGICAL appel1er
+  DATA itapsw /0/
+  DATA appel1er /.TRUE./
+  SAVE itapsw,appel1er
+  !$OMP THREADPRIVATE(appel1er)
+  !$OMP THREADPRIVATE(itapsw)
+  !jq-Introduced for aerosol forcings
+  REAL*8 flag_aer
+  LOGICAL ok_ade, ok_aie    ! use aerosol forcings or not?
+  REAL*8 tauaero(kdlon,kflev,9,2)  ! aerosol optical properties
+  REAL*8 pizaero(kdlon,kflev,9,2)  ! (see aeropt.F)
+  REAL*8 cgaero(kdlon,kflev,9,2)   ! -"-
+  REAL*8 PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
+  REAL*8 POMEGAA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
+  REAL*8 PTOPSWADAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
+  REAL*8 PSOLSWADAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
+  REAL*8 PTOPSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
+  REAL*8 PSOLSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
+  REAL*8 PTOPSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND)
+  REAL*8 PSOLSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND)
+  REAL*8 PTOPSWAERO(KDLON,9)
+  REAL*8 PTOPSW0AERO(KDLON,9)
+  REAL*8 PSOLSWAERO(KDLON,9)
+  REAL*8 PSOLSW0AERO(KDLON,9)
+
+  !jq - Fluxes including aerosol effects
+  REAL*8,ALLOCATABLE,SAVE :: ZFSUPAD_AERO(:,:)
+  !$OMP THREADPRIVATE(ZFSUPAD_AERO)
+  REAL*8,ALLOCATABLE,SAVE :: ZFSDNAD_AERO(:,:)
+  !$OMP THREADPRIVATE(ZFSDNAD_AERO)
+  !jq - Fluxes including aerosol effects
+  REAL*8,ALLOCATABLE,SAVE :: ZFSUPAD0_AERO(:,:)
+  !$OMP THREADPRIVATE(ZFSUPAD0_AERO)
+  REAL*8,ALLOCATABLE,SAVE :: ZFSDNAD0_AERO(:,:)
+  !$OMP THREADPRIVATE(ZFSDNAD0_AERO)
+  REAL*8,ALLOCATABLE,SAVE :: ZFSUPAI_AERO(:,:)
+  !$OMP THREADPRIVATE(ZFSUPAI_AERO)
+  REAL*8,ALLOCATABLE,SAVE :: ZFSDNAI_AERO(:,:)
+  !$OMP THREADPRIVATE(ZFSDNAI_AERO)
+  REAL*8,ALLOCATABLE,SAVE ::  ZFSUP_AERO(:,:,:)
+  !$OMP THREADPRIVATE(ZFSUP_AERO)
+  REAL*8,ALLOCATABLE,SAVE ::  ZFSDN_AERO(:,:,:)
+  !$OMP THREADPRIVATE(ZFSDN_AERO)
+  REAL*8,ALLOCATABLE,SAVE ::  ZFSUP0_AERO(:,:,:)
+  !$OMP THREADPRIVATE(ZFSUP0_AERO)
+  REAL*8,ALLOCATABLE,SAVE ::  ZFSDN0_AERO(:,:,:)
+  !$OMP THREADPRIVATE(ZFSDN0_AERO)
+
+  LOGICAL initialized
+  !rv
+  SAVE flag_aer
+  !$OMP THREADPRIVATE(flag_aer)
+  DATA initialized/.FALSE./
+  SAVE initialized
+  !$OMP THREADPRIVATE(initialized)
+
+
+  IF(.NOT.initialized) THEN
+     flag_aer=0.
+     initialized=.TRUE.
+     ALLOCATE(ZFSUPAD_AERO(KDLON,KFLEV+1))
+     ALLOCATE(ZFSDNAD_AERO(KDLON,KFLEV+1))
+     ALLOCATE(ZFSUPAD0_AERO(KDLON,KFLEV+1))
+     ALLOCATE(ZFSDNAD0_AERO(KDLON,KFLEV+1))
+     ALLOCATE(ZFSUPAI_AERO(KDLON,KFLEV+1))
+     ALLOCATE(ZFSDNAI_AERO(KDLON,KFLEV+1))
+     ALLOCATE(ZFSUP_AERO (KDLON,KFLEV+1,9))
+     ALLOCATE(ZFSDN_AERO (KDLON,KFLEV+1,9))
+     ALLOCATE(ZFSUP0_AERO(KDLON,KFLEV+1,9))
+     ALLOCATE(ZFSDN0_AERO(KDLON,KFLEV+1,9))
+     ZFSUPAD_AERO(:,:)=0.
+     ZFSDNAD_AERO(:,:)=0.
+     ZFSUPAD0_AERO(:,:)=0.
+     ZFSDNAD0_AERO(:,:)=0.
+     ZFSUPAI_AERO(:,:)=0.
+     ZFSDNAI_AERO(:,:)=0.
+     ZFSUP_AERO (:,:,:)=0.
+     ZFSDN_AERO (:,:,:)=0.
+     ZFSUP0_AERO(:,:,:)=0.
+     ZFSDN0_AERO(:,:,:)=0.
+  ENDIF
+  !rv
+
+
+  IF (appel1er) THEN
+     PRINT*, 'SW calling frequency : ', swpas
+     PRINT*, "   In general, it should be 1"
+     appel1er = .FALSE.
+  ENDIF
+  !     ------------------------------------------------------------------
+  IF (MOD(itapsw,swpas).EQ.0) THEN
+
+     DO JK = 1 , KFLEV
+        DO JL = 1, KDLON
+           ZCLDSW0(JL,JK) = 0.0
+           ZOZ(JL,JK) = POZON(JL,JK)*46.6968/RG &
+                *PDP(JL,JK)*(101325.0/PPSOL(JL))
+        ENDDO
+     ENDDO
+
+
+     ! clear-sky:
+     flag_aer=0.0
+     CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+          PRMU0,PFRAC,PTAVE,PWV,&
+          ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+     INU = 1
+     CALL SW1S_LMDAR4(INU,PAER, flag_aer, &
+          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
+          PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
+          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+          ZFD, ZFU)
+     INU = 2
+     CALL SW2S_LMDAR4(INU, PAER, flag_aer, &
+          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
+          ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
+          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+          PWV, PQS,&
+          ZFDOWN, ZFUP)
+     DO JK = 1 , KFLEV+1
+        DO JL = 1, KDLON
+           ZFSUP0(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+           ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ZFSUP0_AERO(JL,JK,1) = ZFSUP0(JL,JK) 
+           ZFSDN0_AERO(JL,JK,1) = ZFSDN0(JL,JK)
+        ENDDO
+     ENDDO
+
+
+     ! cloudy-sky:
+     flag_aer=0.0
+     CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+          PRMU0,PFRAC,PTAVE,PWV,&
+          ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+     INU = 1
+     CALL SW1S_LMDAR4(INU, PAER, flag_aer, &
+          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
+          PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+          ZFD, ZFU)
+     INU = 2
+     CALL SW2S_LMDAR4(INU, PAER, flag_aer, &
+          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
+          ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+          PWV, PQS,&
+          ZFDOWN, ZFUP)
+
+     DO JK = 1 , KFLEV+1
+        DO JL = 1, KDLON
+           ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+           ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ZFSUP_AERO(JL,JK,1) = ZFSUP(JL,JK) 
+           ZFSDN_AERO(JL,JK,1) = ZFSDN(JL,JK)
+        ENDDO
+     ENDDO
+
+     ZFSUP0_AERO(:,:,2) = ZFSUP0_AERO(:,:,1)
+     ZFSDN0_AERO(:,:,2) = ZFSDN0_AERO(:,:,1)
+     ZFSUP_AERO(:,:,2) = ZFSUP_AERO(:,:,1)
+     ZFSDN_AERO(:,:,2) = ZFSDN_AERO(:,:,1)
+
+
+     IF (ok_ade) THEN
+
+        ! clear sky (Anne Cozic 03/07/2007)
+        ! CAS AER (2)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUPAD0_AERO(JL,JK) = ZFSUP0(JL,JK) 
+              ZFSDNAD0_AERO(JL,JK) = ZFSDN0(JL,JK) 
+              ZFSUP0(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+              ZFSUP0_AERO(JL,JK,2) = ZFSUP0(JL,JK) 
+              ZFSDN0_AERO(JL,JK,2) = ZFSDN0(JL,JK)
+           ENDDO
+        ENDDO
+
+        ! cloudy-sky + aerosol dir OB
+        ! ACo AER 
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUPAD_AERO(JL,JK) = ZFSUP(JL,JK) 
+              ZFSDNAD_AERO(JL,JK) = ZFSDN(JL,JK) 
+              ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+              ZFSUP_AERO(JL,JK,2) = ZFSUP(JL,JK) 
+              ZFSDN_AERO(JL,JK,2) = ZFSDN(JL,JK)
+           ENDDO
+        ENDDO
+
+        !CAS NAT
+        ! clear sky
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP0_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! cloudy-sky 
+        ! ACo NAT
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! clear sky (Yves 01/12/2007)
+        ! CAS BC (4)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP0_AERO(JL,JK,4) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! cloudy-sky + aerosol dir OB
+        ! CAS BC (4)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP_AERO(JL,JK,4) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! clear sky (Yves 13/12/2007)
+        ! CAS SO4 (5)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP0_AERO(JL,JK,5) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! cloudy-sky + aerosol dir OB
+        ! CAS SO4 (5)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP_AERO(JL,JK,5) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! clear sky (Yves 13/12/2007)
+        ! CAS POM (6)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP0_AERO(JL,JK,6) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0_AERO(JL,JK,6) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! cloudy-sky + aerosol dir OB
+        ! CAS POM (6)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP_AERO(JL,JK,6) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,6) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! clear sky (Yves 13/12/2007)
+        ! CAS DUST (7)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP0_AERO(JL,JK,7) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0_AERO(JL,JK,7) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! cloudy-sky + aerosol dir OB
+        ! CAS DUST (7)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP_AERO(JL,JK,7) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,7) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! clear sky (Yves 13/12/2007)
+        ! CAS Seasalt SS (8)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP0_AERO(JL,JK,8) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN0_AERO(JL,JK,8) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+        ! cloudy-sky + aerosol dir OB
+        ! CAS Seasalt SS (8)
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUP_AERO(JL,JK,8) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN_AERO(JL,JK,8) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+           ENDDO
+        ENDDO
+
+     ENDIF ! ok_ade
+
+
+     IF (ok_aie) THEN
+        !jq   cloudy-sky + aerosol direct + aerosol indirect
+        flag_aer=1.0
+        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
+             PRMU0,PFRAC,PTAVE,PWV,&
+             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+        INU = 1
+        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
+             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
+             ZFD, ZFU)
+        INU = 2
+        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
+             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
+             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
+             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
+             PWV, PQS,&
+             ZFDOWN, ZFUP)
+        DO JK = 1 , KFLEV+1
+           DO JL = 1, KDLON
+              ZFSUPAI_AERO(JL,JK) = ZFSUP(JL,JK) 
+              ZFSDNAI_AERO(JL,JK) = ZFSDN(JL,JK)          
+              ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
+              ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
+              ZFSUP_AERO(JL,JK,2) = ZFSUP(JL,JK) 
+              ZFSDN_AERO(JL,JK,2) = ZFSDN(JL,JK)
+           ENDDO
+        ENDDO
+     ENDIF ! ok_aie      
+
+     itapsw = 0
+  ENDIF
+  itapsw = itapsw + 1
+
+  DO k = 1, KFLEV
+     kpl1 = k+1
+     DO i = 1, KDLON
+
+        PHEAT(i,k) = -(ZFSUP_AERO(i,kpl1,2)-ZFSUP_AERO(i,k,2)) &
+             -(ZFSDN_AERO(i,k,2)-ZFSDN_AERO(i,kpl1,2))
+        PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
+        PHEAT0(i,k) = -(ZFSUP0_AERO(i,kpl1,2)-ZFSUP0_AERO(i,k,2)) &
+             -(ZFSDN0_AERO(i,k,2)-ZFSDN0_AERO(i,kpl1,2))
+        PHEAT0(i,k) = PHEAT0(i,k) * RDAY*RG/RCPD / PDP(i,k)
+
+     ENDDO
+  ENDDO
+  DO i = 1, KDLON
+     PALBPLA(i) = ZFSUP(i,KFLEV+1)/(ZFSDN(i,KFLEV+1)+1.0e-20)
+     ! clear sky
+     PSOLSW0(i) = ZFSDN0(i,1) - ZFSUP0(i,1)
+     PTOPSW0(i) = ZFSDN0(i,KFLEV+1) - ZFSUP0(i,KFLEV+1)
+
+     PSOLSW(i) = ZFSDN(i,1) - ZFSUP(i,1)
+     PTOPSW(i) = ZFSDN(i,KFLEV+1) - ZFSUP(i,KFLEV+1)
+
+     PSOLSW0AERO(i,:) = ZFSDN0_AERO(i,1,:) - ZFSUP0_AERO(i,1,:)
+     PTOPSW0AERO(i,:) = &
+          ZFSDN0_AERO(i,KFLEV+1,:) - ZFSUP0_AERO(i,KFLEV+1,:)
+
+     PSOLSWAERO(i,:) = ZFSDN_AERO(i,1,:) - ZFSUP_AERO(i,1,:)
+     PTOPSWAERO(i,:) = &
+          ZFSDN_AERO(i,KFLEV+1,:) - ZFSUP_AERO(i,KFLEV+1,:)
+
+     PSOLSWADAERO(i) = ZFSDNAD_AERO(i,1) - ZFSUPAD_AERO(i,1)
+     PTOPSWADAERO(i) = &
+          ZFSDNAD_AERO(i,KFLEV+1) - ZFSUPAD_AERO(i,KFLEV+1)
+
+     PSOLSWAD0AERO(i) = ZFSDNAD0_AERO(i,1) - ZFSUPAD0_AERO(i,1)
+     PTOPSWAD0AERO(i) = &
+          ZFSDNAD0_AERO(i,KFLEV+1) - ZFSUPAD0_AERO(i,KFLEV+1)
+
+     PSOLSWAIAERO(i) = ZFSDNAI_AERO(i,1) - ZFSUPAI_AERO(i,1)
+     PTOPSWAIAERO(i) = &
+          ZFSDNAI_AERO(i,KFLEV+1) - ZFSUPAI_AERO(i,KFLEV+1)
+
+  ENDDO
+END SUBROUTINE SW_AERO
+
