Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 517)
@@ -5,5 +5,7 @@
   subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &
  &                     fact_cldcon, facttemps,ok_newmicro,iflag_cldcon, &
- &                     ratqsbas,ratqshaut,if_ebil)
+ &                     ratqsbas,ratqshaut,if_ebil, &
+ &		       ok_ade, ok_aie, &
+ &                     bl95_b0, bl95_b1)
 
    use IOIPSL
@@ -29,5 +31,7 @@
 ! ok_mensuel: sorties mensuelles
 ! ok_instan:  sorties instantanees
-
+! ok_ade, ok_aie: apply or not aerosol direct and indirect effects
+! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc 
+!
 
 
@@ -35,5 +39,7 @@
   character (len = 6)  :: ocean
   logical              :: ok_veget, ok_newmicro
-  logical              :: ok_journe, ok_mensuel, ok_instan
+  logical              :: ok_journe, ok_mensuel, ok_instan        
+  LOGICAL              :: ok_ade, ok_aie
+  REAL                 :: bl95_b0, bl95_b1
   real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
   integer              :: iflag_cldcon, if_ebil
@@ -89,4 +95,40 @@
   ok_instan = .false.
   call getin('OK_instan', ok_instan)
+!
+!Config Key  = ok_ade
+!Config Desc = Aerosol direct effect or not?
+!Config Def  = .false.
+!Config Help = Used in radlwsw.F
+!
+  ok_ade = .false.
+  call getin('ok_ade', ok_ade)
+
+!
+!Config Key  = ok_aie
+!Config Desc = Aerosol indirect effect or not?
+!Config Def  = .false.
+!Config Help = Used in nuage.F and radlwsw.F
+!
+  ok_aie = .false.
+  call getin('ok_aie', ok_aie)
+
+!
+!Config Key  = bl95_b0
+!Config Desc = Parameter in CDNC-maer link (Boucher&Lohmann 1995)
+!Config Def  = .false.
+!Config Help = Used in nuage.F
+!
+  bl95_b0 = 2.
+  call getin('bl95_b0', bl95_b0)
+
+!Config Key  = bl95_b1
+!Config Desc = Parameter in CDNC-maer link (Boucher&Lohmann 1995)
+!Config Def  = .false.
+!Config Help = Used in nuage.F
+!
+  bl95_b1 = 0.2
+  call getin('bl95_b1', bl95_b1)
+
+!
 !
 !Config Key  = if_ebil
@@ -554,4 +596,8 @@
   write(numout,*)' ksta_ter = ',ksta_ter 
   write(numout,*)' ok_kzmin = ',ok_kzmin 
+  write(numout,*)' ok_ade = ',ok_ade
+  write(numout,*)' ok_aie = ',ok_aie
+  write(numout,*)' bl95_b0 = ',bl95_b0
+  write(numout,*)' bl95_b1 = ',bl95_b1
   write(numout,*)' lev_histhf = ',lev_histhf 
   write(numout,*)' lev_histday = ',lev_histday 
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp_tr.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp_tr.F	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp_tr.F	(revision 517)
@@ -5,5 +5,6 @@
      s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
      s                   frac_impa, frac_nucl,
-     s                   prfl, psfl)
+     s                   prfl, psfl,
+     s                   RHcl) ! relative humidity in clear sky (needed for aer optical properties; aeropt.F)
 
 c
@@ -37,4 +38,8 @@
       REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
       REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
+      
+Cjq   For aerosol opt properties needed (see aeropt.F)
+      REAL RHcl(klon,klev)
+      
 cAA
 c Coeffients de fraction lessivee : pour OFF-LINE
@@ -290,4 +295,11 @@
             rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
             zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
+            
+c--Olivier
+            RHcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+            IF (rneb(i,k) .LE. 0.0) RHcl(i,k)=zq(i)/zqs(i)
+            IF (rneb(i,k) .GE. 1.0) RHcl(i,k)=1.0
+c--fin
+            
          ENDDO
       ELSE
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h	(revision 517)
@@ -277,4 +277,26 @@
      .                "ave(X)", zsto,zout)
 c
+c Effets des aerosols
+c
+c     IF (ok_ade.OR.ok_aie) THEN
+         CALL histdef(nid_mth, "topsad", "ADE at TOA", "W/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "solsad", "ADE at sfc", "W/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "topsai", "AIE at TOA", "W/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "solsai", "AIE at sfc", "W/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "ave(X)", zsto,zout)
+c     endif
+c
+
+c
 c          CALL histdef(nid_mth, "frtu", "Zonal wind stress", "Pa",
 c    .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
@@ -693,5 +715,29 @@
      .                "ave(X)", zsto,zout)
          ENDIF
-C
+c
+c Effets des aerosols
+c
+c     IF (ok_ade.OR.ok_aie) THEN
+         CALL histdef(nid_mth, "re", "CDR", "um",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "redenom", "CDR denominator", "-",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "tau", "cloud opt thickness", "-",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "taupi", "cloud opt thickn. (pi)", "-",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+c     endif
+c
+         CALL histdef(nid_mth, "ozone", "Ozone concentration", "-",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+c
          if (nqmax.GE.3) THEN
          DO iq=1,nqmax-2
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/newmicro.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/newmicro.F	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/newmicro.F	(revision 517)
@@ -1,8 +1,10 @@
       SUBROUTINE newmicro (paprs, pplay,ok_newmicro,
      .                  t, pqlwp, pclc, pcltau, pclemi,
-cIM    .                  pch, pcl, pcm, pct, pctlwp)
      .                  pch, pcl, pcm, pct, pctlwp,
-     .                  xflwp, xfiwp, xflwc, xfiwc)
-
+     s                  xflwp, xfiwp, xflwc, xfiwc,
+     e                  ok_aie, 
+     e                  sulfate, sulfate_pi, 
+     e                  bl95_b0, bl95_b1,
+     s                  cldtaupi, re, fl)
       IMPLICIT none
 c======================================================================
@@ -15,4 +17,17 @@
 c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
 c 
+c ok_aie--input-L-apply aerosol indirect effect or not
+c sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
+c sulfate_pi-input-R-dito, pre-industrial value
+c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
+c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
+c      
+c cldtaupi-output-R-pre-industrial value of cloud optical thickness, 
+c                   needed for the diagnostics of the aerosol indirect 
+c                   radiative forcing (see radlwsw)
+c re------output-R-Cloud droplet effective radius multiplied by fl [um]
+c fl------output-R-Denominator to re, introduced to avoid problems in
+c                  the averaging of the output. fl is the fraction of liquid
+c                  water clouds within a grid cell           
 c pcltau--output-R-epaisseur optique des nuages
 c pclemi--output-R-emissivite des nuages (0 a 1)
@@ -66,5 +81,23 @@
       parameter (DF=1.66) ! diffusivity factor
 c sb --
-
+cjq for the aerosol indirect effect
+cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
+cjq      
+      LOGICAL ok_aie            ! Apply AIE or not?
+      LOGICAL ok_a1lwpdep       ! a1 LWP dependent?
+      
+      REAL sulfate(klon, klev)  ! sulfate aerosol mass concentration [ug m-3]
+      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
+      REAL re(klon, klev)       ! cloud droplet effective radius [um]
+      REAL sulfate_pi(klon, klev)  ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
+      REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
+      REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
+      
+      REAL fl(klon, klev)       ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
+      
+      REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
+      
+      REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
+cjq-end    
 c
 c Calculer l'epaisseur optique et l'emmissivite des nuages
@@ -119,6 +152,68 @@
 
 c for liquid water clouds: 
+         IF (ok_aie) THEN
+            ! Formula "D" of Boucher and Lohmann, Tellus, 1995
+            !             
+            cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
+     .           log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
+            ! Cloud droplet number concentration (CDNC) is restricted
+            ! to be within [20, 1000 cm^3]
+            ! 
+            cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
+            !
+            !
+            cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
+     .           log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
+            cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
+            !            
+            !
+            ! air density: pplay(i,k) / (RD * zT(i,k)) 
+            ! factor 1.1: derive effective radius from volume-mean radius
+            ! factor 1000 is the water density
+            ! _chaud means that this is the CDR for liquid water clouds
+            !
+            rad_chaud = 
+     .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) )  
+     .               / (4./3. * RPI * 1000. * cdnc(i,k)) )**(1./3.)
+            !
+            ! Convert to um. CDR shall be at least 3 um.
+            !
+c           rad_chaud = MAX(rad_chaud*1.e6, 3.) 
+            rad_chaud = MAX(rad_chaud*1.e6, 5.) 
+            
+            ! Pre-industrial cloud opt thickness
+            !
+            ! "radius" is calculated as rad_chaud above (plus the 
+            ! ice cloud contribution) but using cdnc_pi instead of
+            ! cdnc.
+            radius = 
+     .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) )  
+     .               / (4./3. * RPI * 1000. * cdnc_pi(i,k)) )**(1./3.)
+            radius = MAX(radius*1.e6, 3.) 
+            
+            tc = t(i,k)-273.15
+            rei = 0.71*tc + 61.29 
+            if (tc.le.-81.4) rei = 3.5 
+            if (zflwp(i).eq.0.) radius = 1. 
+            if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1. 
+            cldtaupi(i,k) = 3.0/2.0 * zflwp(i) / radius
+     .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
+         ENDIF                  ! ok_aie
+         ! For output diagnostics
+         !
+         ! Cloud droplet effective radius [um]
+         !
+         ! we multiply here with f * xl (fraction of liquid water
+         ! clouds in the grid cell) to avoid problems in the
+         ! averaging of the output.
+         ! In the output of IOIPSL, derive the real cloud droplet 
+         ! effective radius as re/fl
+         !
+         fl(i,k) = pclc(i,k)*(1.-zfice)            
+         re(i,k) = rad_chaud*fl(i,k)
+            
+c-jq end         
+         
          rel = rad_chaud
-
 c for ice clouds: as a function of the ambiant temperature
 c [formula used by Iacobellis and Somerville (2000), with an 
@@ -156,4 +251,7 @@
          IF (lo) pcltau(i,k) = 0.0
          IF (lo) pclemi(i,k) = 0.0
+         
+         IF (lo) cldtaupi(i,k) = 0.0
+         IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)            
       ENDDO
       ENDDO
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/nuage.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/nuage.F	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/nuage.F	(revision 517)
@@ -1,5 +1,9 @@
       SUBROUTINE nuage (paprs, pplay,
      .                  t, pqlwp, pclc, pcltau, pclemi,
-     .                  pch, pcl, pcm, pct, pctlwp)
+     .                  pch, pcl, pcm, pct, pctlwp,
+     e                  ok_aie,
+     e                  sulfate, sulfate_pi, 
+     e                  bl95_b0, bl95_b1,
+     s                  cldtaupi, re, fl)
       IMPLICIT none
 c======================================================================
@@ -11,4 +15,17 @@
 c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
 c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
+c ok_aie--input-L-apply aerosol indirect effect or not
+c sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
+c sulfate_pi-input-R-dito, pre-industrial value
+c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
+c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
+c      
+c cldtaupi-output-R-pre-industrial value of cloud optical thickness, 
+c                   needed for the diagnostics of the aerosol indirect 
+c                   radiative forcing (see radlwsw)
+c re------output-R-Cloud droplet effective radius multiplied by fl [um]
+c fl------output-R-Denominator to re, introduced to avoid problems in
+c                  the averaging of the output. fl is the fraction of liquid
+c                  water clouds within a grid cell      
 c 
 c pcltau--output-R-epaisseur optique des nuages
@@ -20,5 +37,4 @@
 #include "dimensions.h"
 #include "dimphy.h"
-#include "nuage.h"
       REAL paprs(klon,klev+1), pplay(klon,klev)
       REAL t(klon,klev)
@@ -38,6 +54,6 @@
       REAL zflwp, zradef, zfice, zmsac
 c
-      REAL radius, rad_chaud
-ccc      PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
+      REAL radius, rad_froid, rad_chaud, rad_chau1, rad_chau2
+      PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
 ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
 c sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
@@ -48,4 +64,24 @@
       INTEGER nexpo ! exponentiel pour glace/eau
       PARAMETER (nexpo=6)
+      
+cjq for the aerosol indirect effect
+cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
+cjq      
+      LOGICAL ok_aie            ! Apply AIE or not?
+      
+      REAL sulfate(klon, klev)  ! sulfate aerosol mass concentration [ug m-3]
+      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
+      REAL re(klon, klev)       ! cloud droplet effective radius [um]
+      REAL sulfate_pi(klon, klev)  ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
+      REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
+      REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
+      
+      REAL fl(klon, klev)  ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
+      
+      REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
+      
+      REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
+cjq-end      
+      
 ccc      PARAMETER (nexpo=1)
 c
@@ -56,4 +92,5 @@
          rad_chaud = rad_chau1
          IF (k.LE.3) rad_chaud = rad_chau2
+            
          pclc(i,k) = MAX(pclc(i,k), seuil_neb)
          zflwp = 1000.*pqlwp(i,k)/RG/pclc(i,k)
@@ -62,4 +99,57 @@
          zfice = MIN(MAX(zfice,0.0),1.0)
          zfice = zfice**nexpo
+         
+         IF (ok_aie) THEN
+            ! Formula "D" of Boucher and Lohmann, Tellus, 1995
+            !             
+            cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
+     .           log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
+            ! Cloud droplet number concentration (CDNC) is restricted
+            ! to be within [20, 1000 cm^3]
+            ! 
+            cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
+            cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
+     .           log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
+            cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
+            !            
+            !
+            ! air density: pplay(i,k) / (RD * zT(i,k)) 
+            ! factor 1.1: derive effective radius from volume-mean radius
+            ! factor 1000 is the water density
+            ! _chaud means that this is the CDR for liquid water clouds
+            !
+            rad_chaud = 
+     .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) )  
+     .               / (4./3. * RPI * 1000. * cdnc(i,k)) )**(1./3.)
+            !
+            ! Convert to um. CDR shall be at least 3 um.
+            !
+            rad_chaud = MAX(rad_chaud*1.e6, 3.) 
+            
+            ! For output diagnostics
+            !
+            ! Cloud droplet effective radius [um]
+            !
+            ! we multiply here with f * xl (fraction of liquid water
+            ! clouds in the grid cell) to avoid problems in the
+            ! averaging of the output.
+            ! In the output of IOIPSL, derive the real cloud droplet 
+            ! effective radius as re/fl
+            !
+            fl(i,k) = pclc(i,k)*(1.-zfice)            
+            re(i,k) = rad_chaud*fl(i,k)
+            
+            ! Pre-industrial cloud opt thickness
+            !
+            ! "radius" is calculated as rad_chaud above (plus the 
+            ! ice cloud contribution) but using cdnc_pi instead of
+            ! cdnc.
+            radius = MAX(1.1e6 * ( (pqlwp(i,k)*pplay(i,k)/(RD*T(i,k)))  
+     .                / (4./3.*RPI*1000.*cdnc_pi(i,k)) )**(1./3.), 
+     .               3.) * (1.-zfice) + rad_froid * zfice           
+            cldtaupi(i,k) = 3.0/2.0 * zflwp / radius
+     .           
+         ENDIF                  ! ok_aie
+         
          radius = rad_chaud * (1.-zfice) + rad_froid * zfice
          coef = coef_chau * (1.-zfice) + coef_froi * zfice
@@ -70,4 +160,6 @@
          IF (lo) pcltau(i,k) = 0.0
          IF (lo) pclemi(i,k) = 0.0
+         
+         IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)            
       ENDDO
       ENDDO
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 517)
@@ -1,3 +1,3 @@
-c
+C
 c $Header$
 c
@@ -143,4 +143,8 @@
       LOGICAL ok_region ! sortir le fichier regional
       PARAMETER (ok_region=.FALSE.)
+c
+c
+      LOGICAL ok_polder ! sortir échantillonné de manière POLDER
+      save ok_polder
 c======================================================================
 c
@@ -877,6 +881,6 @@
       REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
 c
-      INTEGER nid_day, nid_mth, nid_ins, nid_nmc
-      SAVE nid_day, nid_mth, nid_ins, nid_nmc
+      INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_pol
+      SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_pol
 c
       INTEGER nhori, nvert
@@ -928,4 +932,38 @@
       REAL zu10m(klon), zv10m(klon)           !vents a 10m moyennes s/1 maille
       CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
+cjq   Aerosol effects (Johannes Quaas, 27/11/2003)
+      REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
+      REAL sulfate_pi(klon, klev) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
+      SAVE sulfate_pi
+
+      REAL cldtaupi(klon,klev)  ! Cloud optical thickness for pre-industrial (pi) aerosols
+
+      REAL re(klon, klev)       ! Cloud droplet effective radius
+      REAL fl(klon, klev)  ! denominator of re
+
+      REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
+
+      ! Aerosol optical properties
+      REAL tau_ae(klon,klev,2), piz_ae(klon,klev,2)
+      REAL cg_ae(klon,klev,2)
+
+      REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
+      ! ok_ade=T -ADE=topswad-topsw
+
+      REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
+      ! ok_aie=T ->
+      !        ok_ade=T -AIE=topswai-topswad
+      !        ok_ade=F -AIE=topswai-topsw
+
+      ! For POLDER swath
+      INTEGER pyr, pmo, pday    ! Year, month and day
+      INTEGER poldermask(klon)  ! POLDER swath mask (0 or 1)
+
+      REAL aerindex(klon)       ! POLDER aerosol index
+     
+      ! Parameters
+      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
+      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
+cjq-end
 c
 c Declaration des constantes et des fonctions thermodynamiques
@@ -991,5 +1029,7 @@
          call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
      .                  ok_instan, fact_cldcon, facttemps,ok_newmicro,
-     .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil)
+     .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
+     .                  ok_ade, ok_aie, 
+     .                  bl95_b0, bl95_b1)
 cIM  .                  , RI0)
 
@@ -1145,5 +1185,5 @@
 
 #undef histmthNMC
-#define histmthNMC
+cccccccc#define histmthNMC
 #ifdef histmthNMC
 #include "ini_histmthNMC.h"
@@ -1151,4 +1191,5 @@
 
 #include "ini_histins.h"
+#include "ini_histpol.h"
 
 #ifdef histREGDYN
@@ -1912,5 +1953,9 @@
      .            t_seri, convliq, convfra, dtau_c, dem_c,
      .            cldh_c, cldl_c, cldm_c, cldt_c, cldq_c,
-     .            flwp_c, fiwp_c, flwc_c, fiwc_c)
+     .            flwp_c, fiwp_c, flwc_c, fiwc_c,
+     e            ok_aie,
+     e            sulfate, sulfate_pi,
+     e            bl95_b0, bl95_b1,
+     s            cldtaupi, re, fl)
 c
 cIM calcul tau. emi nuages startiformes
@@ -1918,5 +1963,9 @@
      .            t_seri, cldliq, cldfra, dtau_s, dem_s,
      .            cldh_s, cldl_s, cldm_s, cldt_s, cldq_s,
-     .            flwp_s, fiwp_s, flwc_s, fiwc_s)
+     .            flwp_s, fiwp_s, flwc_s, fiwc_s,
+     e            ok_aie,
+     e            sulfate, sulfate_pi,
+     e            bl95_b0, bl95_b1,
+     s            cldtaupi, re, fl)
 c
       cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
@@ -2185,5 +2234,17 @@
       ENDDO
       ENDDO
-c
+cjq - introduce the aerosol direct and first indirect radiative forcings
+cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
+      IF (ok_ade.OR.ok_aie) 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
+
+c     
 c Calculer les parametres optiques des nuages et quelques
 c parametres pour diagnostiques:
@@ -2193,9 +2254,18 @@
      .            t_seri, cldliq, cldfra, cldtau, cldemi,
      .            cldh, cldl, cldm, cldt, cldq,
-     .            flwp, fiwp, flwc, fiwc)
+     .            flwp, fiwp, flwc, fiwc,
+     e            ok_aie,
+     e            sulfate, sulfate_pi,
+     e            bl95_b0, bl95_b1,
+     s            cldtaupi, re, fl)
       else
       CALL nuage (paprs, pplay,
      .            t_seri, cldliq, cldfra, cldtau, cldemi,
-     .            cldh, cldl, cldm, cldt, cldq)
+     .            cldh, cldl, cldm, cldt, cldq,
+     e            ok_aie,
+     e            sulfate, sulfate_pi,
+     e            bl95_b0, bl95_b1,
+     s            cldtaupi, re, fl)
+      
       endif
 c
@@ -2233,5 +2303,10 @@
      s             topsw0,toplw0,solsw0,sollw0,
      s             lwdn0, lwdn, lwup0, lwup, 
-     s             swdn0, swdn, swup0, swup     )
+     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) ! ="=
       itaprad = 0
       ENDIF
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 517)
@@ -1,16 +1,16 @@
-cIM   SUBROUTINE radlwsw(dist, rmu0, fract, co2_ppm, solaire,
       SUBROUTINE radlwsw(dist, rmu0, fract, 
      .                  paprs, pplay,tsol,albedo, alblw, t,q,wo,
-     .                  cldfra, cldemi, cldtau,
+     .                  cldfra, cldemi, cldtaupd,
      .                  heat,heat0,cool,cool0,radsol,albpla,
      .                  topsw,toplw,solsw,sollw,
      .                  sollwdown,
-cIM  .                  sollwdown, sollwdownclr,
-cIM  .                  toplwdown, toplwdownclr,
      .                  topsw0,toplw0,solsw0,sollw0,
-cIM BEG
      .                  lwdn0, lwdn, lwup0, lwup,
-cIM END
-     .                  swdn0, swdn, swup0, swup    )
+     .                  swdn0, swdn, swup0, swup,
+     .                  ok_ade, ok_aie,
+     .                  tau_ae, piz_ae, cg_ae,
+     .                  topswad, solswad,
+     .                  cldtaupi, topswai, solswai)
+c      
       IMPLICIT none
 c======================================================================
@@ -31,6 +31,13 @@
 c wo-------input-R- contenu en ozone (en cm.atm)
 c cldfra---input-R- fraction nuageuse (entre 0 et 1)
-c cldtau---input-R- epaisseur optique des nuages dans le visible
+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)
@@ -42,4 +49,20 @@
 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======================================================================
 #include "dimensions.h"
@@ -56,5 +79,5 @@
       real albedo(klon), alblw(klon), tsol(klon)
       real t(klon,klev), q(klon,klev), wo(klon,klev)
-      real cldfra(klon,klev), cldemi(klon,klev), cldtau(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)
@@ -123,5 +146,24 @@
       REAL lwup(klon,kflev+1),lwup0(klon,kflev+1)
 cIM END
-c---------------------------------------------------------------
+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
+      
+c
+c-------------------------------------------
       nb_gr = klon / kdlon
       IF (nb_gr*kdlon .NE. klon) THEN
@@ -202,10 +244,21 @@
          PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
          PCLDSW(i,k) = cldfra(iof+i,k)
-         PTAU(i,1,k) = MAX(cldtau(iof+i,k), 1.0e-05)! 1e-12 serait instable
-         PTAU(i,2,k) = MAX(cldtau(iof+i,k), 1.0e-05)! pour 32-bit machines
+         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
@@ -222,4 +275,15 @@
          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
@@ -247,5 +311,10 @@
      S        zheat, zheat0,
      S        zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,
-     S        ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
+     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?
+
 c======================================================================
       DO i = 1, kdlon
@@ -292,4 +361,29 @@
 c        swup  ( iof+i,2)   = ZFSUP  ( i,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
@@ -321,5 +415,10 @@
      S              PHEAT, PHEAT0,
      S              PALBPLA,PTOPSW,PSOLSW,PTOPSW0,PSOLSW0,
-     S              ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
+     S              ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,
+     S              tauae, pizae, cgae,
+     s              PTAUA, POMEGAA,
+     S              PTOPSWAD,PSOLSWAD,PTOPSWAI,PSOLSWAI,
+     J              ok_ade, ok_aie )
+      
       IMPLICIT none
 
@@ -358,4 +457,5 @@
 C        ORIGINAL : 89-07-14
 C        95-01-01   J.-J. MORCRETTE  Direct/Diffuse Albedo
+c        03-11-27   J. QUAAS Introduce aerosol forcings (based on BOUCHER)
 C     ------------------------------------------------------------------
 C
@@ -426,4 +526,24 @@
       DATA itapsw /0/
       DATA appel1er /.TRUE./
+cjq-Introduced for aerosol forcings
+      real*8 flag_aer
+      logical ok_ade, ok_aie    ! use aerosol forcings or not?
+      real*8 tauae(kdlon,kflev,2)  ! aerosol optical properties
+      real*8 pizae(kdlon,kflev,2)  ! (see aeropt.F)
+      real*8 cgae(kdlon,kflev,2)   ! -"-
+      REAL*8 PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
+      REAL*8 POMEGAA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
+      REAL*8 PTOPSWAD(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
+      REAL*8 PSOLSWAD(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
+      REAL*8 PTOPSWAI(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND)
+      REAL*8 PSOLSWAI(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND)
+cjq - Fluxes including aerosol effects
+      REAL*8 ZFSUPAD(KDLON,KFLEV+1)
+      REAL*8 ZFSDNAD(KDLON,KFLEV+1)
+      REAL*8 ZFSUPAI(KDLON,KFLEV+1)
+      REAL*8 ZFSDNAI(KDLON,KFLEV+1)
+      SAVE ZFSUPAD, ZFSDNAD, ZFSUPAI, ZFSDNAI ! aerosol fluxes
+cjq-end
+      
 c
       IF (appel1er) THEN
@@ -451,10 +571,12 @@
       INU = 1
       CALL SW1S(INU,
-     S     PAER, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,
      S     ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,
      S     ZFD, ZFU)
       INU = 2
       CALL SW2S(INU,
-     S     PAER, ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,
      S     ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,
      S     PWV, PQS,
@@ -466,6 +588,6 @@
       ENDDO
       ENDDO
-c cloudy-sky:
-cIM ctes ds clesphys.h   CALL SWU(PSCT,RCO2,PCLDSW,PPMB,PPSOL,
+      
+      flag_aer=0.0
       CALL SWU(PSCT,PCLDSW,PPMB,PPSOL,
      S         PRMU0,PFRAC,PTAVE,PWV,
@@ -473,13 +595,18 @@
       INU = 1
       CALL SW1S(INU,
-     S     PAER, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
      S     ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,
      S     ZFD, ZFU)
       INU = 2
       CALL SW2S(INU,
-     S     PAER, ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
      S     ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,
      S     PWV, PQS,
      S    ZFDOWN, ZFUP)
+
+c cloudy-sky:
+      
       DO JK = 1 , KFLEV+1
       DO JL = 1, KDLON
@@ -488,5 +615,68 @@
       ENDDO
       ENDDO
+      
+c      
+      IF (ok_ade) THEN
 c
+c cloudy-sky + aerosol dir OB
+      flag_aer=1.0
+      CALL SWU(PSCT,PCLDSW,PPMB,PPSOL,
+     S         PRMU0,PFRAC,PTAVE,PWV,
+     S         ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+      INU = 1
+      CALL SW1S(INU,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
+     S     ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,
+     S     ZFD, ZFU)
+      INU = 2
+      CALL SW2S(INU,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
+     S     ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,
+     S     PWV, PQS,
+     S    ZFDOWN, ZFUP)
+      DO JK = 1 , KFLEV+1
+      DO JL = 1, KDLON
+         ZFSUPAD(JL,JK) = ZFSUP(JL,JK) 
+         ZFSDNAD(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)
+      ENDDO
+      ENDDO 
+      
+      ENDIF ! ok_ade
+      
+      IF (ok_aie) THEN
+         
+cjq   cloudy-sky + aerosol direct + aerosol indirect
+      flag_aer=1.0
+      CALL SWU(PSCT,PCLDSW,PPMB,PPSOL,
+     S         PRMU0,PFRAC,PTAVE,PWV,
+     S         ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
+      INU = 1
+      CALL SW1S(INU,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
+     S     ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,
+     S     ZFD, ZFU)
+      INU = 2
+      CALL SW2S(INU,
+     S     PAER, flag_aer, tauae, pizae, cgae,
+     S     ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,
+     S     ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,
+     S     PWV, PQS,
+     S    ZFDOWN, ZFUP)
+      DO JK = 1 , KFLEV+1
+      DO JL = 1, KDLON
+         ZFSUPAI(JL,JK) = ZFSUP(JL,JK) 
+         ZFSDNAI(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)
+      ENDDO
+      ENDDO
+      ENDIF ! ok_aie      
+cjq -end
+      
       itapsw = 0
       ENDIF
@@ -512,4 +702,11 @@
          PSOLSW0(i) = ZFSDN0(i,1) - ZFSUP0(i,1)
          PTOPSW0(i) = ZFSDN0(i,KFLEV+1) - ZFSUP0(i,KFLEV+1)
+c-OB
+         PSOLSWAD(i) = ZFSDNAD(i,1) - ZFSUPAD(i,1)
+         PTOPSWAD(i) = ZFSDNAD(i,KFLEV+1) - ZFSUPAD(i,KFLEV+1)
+c
+         PSOLSWAI(i) = ZFSDNAI(i,1) - ZFSUPAI(i,1)
+         PTOPSWAI(i) = ZFSDNAI(i,KFLEV+1) - ZFSUPAI(i,KFLEV+1)
+c-fin 
       ENDDO
 C
@@ -707,5 +904,6 @@
       END
       SUBROUTINE SW1S ( KNU
-     S  ,  PAER  , PALBD , PALBP, PCG  , PCLD , PCLEAR, PCLDSW
+     S  ,  PAER  , flag_aer, tauae, pizae, cgae
+     S  ,  PALBD , PALBP, PCG  , PCLD , PCLEAR, PCLDSW
      S  ,  PDSIG , POMEGA, POZ  , PRMU , PSEC , PTAU  , PUD  
      S  ,  PFD   , PFU)
@@ -748,4 +946,9 @@
 C
       INTEGER KNU
+c-OB
+      real*8 flag_aer
+      real*8 tauae(kdlon,kflev,2)
+      real*8 pizae(kdlon,kflev,2)
+      real*8 cgae(kdlon,kflev,2)
       REAL*8 PAER(KDLON,KFLEV,5)
       REAL*8 PALBD(KDLON,2)
@@ -839,5 +1042,6 @@
 C
       CALL SWCLR ( KNU
-     S  , PAER   , PALBP , PDSIG , ZRAYL, PSEC
+     S  , PAER   , flag_aer, tauae, pizae, cgae
+     S  , PALBP  , PDSIG , ZRAYL, PSEC
      S  , ZCGAZ  , ZPIZAZ, ZRAY1 , ZRAY2, ZREFZ, ZRJ0
      S  , ZRK0   , ZRMU0 , ZTAUAZ, ZTRA1, ZTRA2)
@@ -939,5 +1143,6 @@
       END
       SUBROUTINE SW2S ( KNU
-     S  ,  PAER  ,PAKI, PALBD, PALBP, PCG   , PCLD, PCLEAR, PCLDSW
+     S  ,  PAER  , flag_aer, tauae, pizae, cgae
+     S  ,  PAKI, PALBD, PALBP, PCG   , PCLD, PCLEAR, PCLDSW
      S  ,  PDSIG ,POMEGA,POZ , PRMU , PSEC  , PTAU
      S  ,  PUD   ,PWV , PQS
@@ -986,4 +1191,9 @@
 C
       INTEGER KNU
+c-OB
+      real*8 flag_aer
+      real*8 tauae(kdlon,kflev,2)
+      real*8 pizae(kdlon,kflev,2)
+      real*8 cgae(kdlon,kflev,2)
       REAL*8 PAER(KDLON,KFLEV,5)
       REAL*8 PAKI(KDLON,2)
@@ -1107,5 +1317,6 @@
 C
       CALL SWCLR ( KNU
-     S  , PAER   , PALBP , PDSIG , ZRAYL, PSEC 
+     S  , PAER   , flag_aer, tauae, pizae, cgae
+     S  , PALBP  , PDSIG , ZRAYL, PSEC 
      S  , ZCGAZ  , ZPIZAZ, ZRAY1 , ZRAY2, ZREFZ, ZRJ0
      S  , ZRK0   , ZRMU0 , ZTAUAZ, ZTRA1, ZTRA2)
@@ -1479,5 +1690,6 @@
       END
       SUBROUTINE SWCLR  ( KNU
-     S  , PAER  , PALBP , PDSIG , PRAYL , PSEC
+     S  , PAER  , flag_aer, tauae, pizae, cgae
+     S  , PALBP , PDSIG , PRAYL , PSEC
      S  , PCGAZ , PPIZAZ, PRAY1 , PRAY2 , PREFZ , PRJ  
      S  , PRK   , PRMU0 , PTAUAZ, PTRA1 , PTRA2                   )
@@ -1512,4 +1724,9 @@
 C
       INTEGER KNU
+c-OB
+      real*8 flag_aer
+      real*8 tauae(kdlon,kflev,2)
+      real*8 pizae(kdlon,kflev,2)
+      real*8 cgae(kdlon,kflev,2)
       REAL*8 PAER(KDLON,KFLEV,5)
       REAL*8 PALBP(KDLON,2)
@@ -1576,24 +1793,33 @@
 C
       DO 108 JK = 1 , KFLEV
-      DO 104 JL = 1, KDLON
-      PCGAZ(JL,JK) = 0.
-      PPIZAZ(JL,JK) =  0.
-      PTAUAZ(JL,JK) = 0.
- 104  CONTINUE
-      DO 106 JAE=1,5
+c-OB
+c      DO 104 JL = 1, KDLON
+c      PCGAZ(JL,JK) = 0.
+c      PPIZAZ(JL,JK) =  0.
+c      PTAUAZ(JL,JK) = 0.
+c 104  CONTINUE
+c-OB
+c      DO 106 JAE=1,5
+c      DO 105 JL = 1, KDLON
+c      PTAUAZ(JL,JK)=PTAUAZ(JL,JK)
+c     S        +PAER(JL,JK,JAE)*TAUA(KNU,JAE)
+c      PPIZAZ(JL,JK)=PPIZAZ(JL,JK)+PAER(JL,JK,JAE)
+c     S        * TAUA(KNU,JAE)*RPIZA(KNU,JAE)
+c      PCGAZ(JL,JK) =  PCGAZ(JL,JK) +PAER(JL,JK,JAE)
+c     S        * TAUA(KNU,JAE)*RPIZA(KNU,JAE)*RCGA(KNU,JAE)
+c 105  CONTINUE
+c 106  CONTINUE
+c-OB
       DO 105 JL = 1, KDLON
-      PTAUAZ(JL,JK)=PTAUAZ(JL,JK)
-     S        +PAER(JL,JK,JAE)*TAUA(KNU,JAE)
-      PPIZAZ(JL,JK)=PPIZAZ(JL,JK)+PAER(JL,JK,JAE)
-     S        * TAUA(KNU,JAE)*RPIZA(KNU,JAE)
-      PCGAZ(JL,JK) =  PCGAZ(JL,JK) +PAER(JL,JK,JAE)
-     S        * TAUA(KNU,JAE)*RPIZA(KNU,JAE)*RCGA(KNU,JAE)
+      PTAUAZ(JL,JK)=flag_aer * tauae(JL,JK,KNU)
+      PPIZAZ(JL,JK)=flag_aer * pizae(JL,JK,KNU)
+      PCGAZ (JL,JK)=flag_aer * cgae(JL,JK,KNU)
  105  CONTINUE
- 106  CONTINUE
-C
+C
+      IF (flag_aer.GT.0) THEN
+c-OB
       DO 107 JL = 1, KDLON
-      IF (KAER.NE.0) THEN
-         PCGAZ(JL,JK)=PCGAZ(JL,JK)/PPIZAZ(JL,JK)
-         PPIZAZ(JL,JK)=PPIZAZ(JL,JK)/PTAUAZ(JL,JK)
+c         PCGAZ(JL,JK)=PCGAZ(JL,JK)/PPIZAZ(JL,JK)
+c         PPIZAZ(JL,JK)=PPIZAZ(JL,JK)/PTAUAZ(JL,JK)
          ZTRAY = PRAYL(JL) * PDSIG(JL,JK)
          ZRATIO = ZTRAY / (ZTRAY + PTAUAZ(JL,JK))
@@ -1604,11 +1830,14 @@
          PPIZAZ(JL,JK) =ZRATIO+(1.-ZRATIO)*PPIZAZ(JL,JK)*(1.-ZFF)
      S                       / (1. - PPIZAZ(JL,JK) * ZFF)
+ 107  CONTINUE
       ELSE
+      DO JL = 1, KDLON
          ZTRAY = PRAYL(JL) * PDSIG(JL,JK)
          PTAUAZ(JL,JK) = ZTRAY
          PCGAZ(JL,JK) = 0.
          PPIZAZ(JL,JK) = 1.-REPSCT
-      END IF
- 107  CONTINUE
+      END DO
+      END IF   ! check flag_aer
+c     107  CONTINUE
 c      PRINT 9107,JK,((PAER(JL,JK,JAE),JAE=1,5)
 c     $ ,PTAUAZ(JL,JK),PPIZAZ(JL,JK),PCGAZ(JL,JK),JL=1,KDLON)
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h	(revision 515)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h	(revision 517)
@@ -262,4 +262,29 @@
      $               iim*jjmp1,ndex2d)
 c
+c
+c effets des aerosols
+c
+c     IF (ok_ade.OR.ok_aie) THEN
+      zx_tmp_fi2d(1:klon) = topswai(1:klon) - topswad(1:klon)
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topswad,zx_tmp_2d)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1,zx_tmp_fi2d ,zx_tmp_2d)
+      CALL histwrite(nid_mth,"topsad",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      zx_tmp_fi2d(1:klon) = solswai(1:klon) - solswad(1:klon)
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solswad,zx_tmp_2d)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1,zx_tmp_fi2d ,zx_tmp_2d)
+      CALL histwrite(nid_mth,"solsad",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      zx_tmp_fi2d(1:klon) = topsw(1:klon) - topswai(1:klon)
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topswai,zx_tmp_2d)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1,zx_tmp_fi2d ,zx_tmp_2d)
+      CALL histwrite(nid_mth,"topsai",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      zx_tmp_fi2d(1:klon) = solsw(1:klon) - solswai(1:klon)
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solswai,zx_tmp_2d)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1,zx_tmp_fi2d ,zx_tmp_2d)
+      CALL histwrite(nid_mth,"solsai",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c     endif
+c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
       CALL histwrite(nid_mth,"bils",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
@@ -712,4 +737,25 @@
       ENDIF
 C
+c
+c effets des aerosols
+c
+c     IF (ok_ade.OR.ok_aie) THEN
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, re, zx_tmp_3d)
+      CALL histwrite(nid_mth,"re",itau_w,zx_tmp_3d,
+     .                                   iim*jjmp1*klev,ndex3d)
+c
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, fl, zx_tmp_3d)
+      CALL histwrite(nid_mth,"redenom",itau_w,zx_tmp_3d,
+     .                                   iim*jjmp1*klev,ndex3d)
+c
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldtau, zx_tmp_3d)
+      CALL histwrite(nid_mth,"tau",itau_w,zx_tmp_3d,
+     .                                   iim*jjmp1*klev,ndex3d)
+c
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldtaupi, zx_tmp_3d)
+      CALL histwrite(nid_mth,"taupi",itau_w,zx_tmp_3d,
+     .                                   iim*jjmp1*klev,ndex3d)
+c     endif
+c
       IF (nqmax.GE.3) THEN
       DO iq=1,nqmax-2
