Index: trunk/LMDZ.PLUTO.old/libf/phypluto/aerosol_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO.old/libf/phypluto/aerosol_mod.F90	(revision 3357)
+++ trunk/LMDZ.PLUTO.old/libf/phypluto/aerosol_mod.F90	(revision 3373)
@@ -55,4 +55,5 @@
       parameter(Nfine=701)
       character(len=100) :: file_path
+      character(len=100) :: file_name
       real,save :: levdat(Nfine),densdat(Nfine)
 
Index: trunk/LMDZ.PLUTO.old/libf/phypluto/callcorrk.F
===================================================================
--- trunk/LMDZ.PLUTO.old/libf/phypluto/callcorrk.F	(revision 3357)
+++ trunk/LMDZ.PLUTO.old/libf/phypluto/callcorrk.F	(revision 3373)
@@ -317,6 +317,6 @@
 !     Prepare NON LTE correction in Pluto atmosphere
       IF (nlte) then
-    !     CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
-    !  &             eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
+        CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
+     &             eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
       ENDIF
 c     Net atmospheric radiative cooling rate from C2H2 (K.s-1):
Index: trunk/LMDZ.PLUTO.old/libf/phypluto/inifis.F
===================================================================
--- trunk/LMDZ.PLUTO.old/libf/phypluto/inifis.F	(revision 3357)
+++ trunk/LMDZ.PLUTO.old/libf/phypluto/inifis.F	(revision 3373)
@@ -135,13 +135,13 @@
          hazerad_file="hazerad.txt"  ! default file
          call getin("hazerad_file",hazerad_file)
-         write(*,*) trim(rname)//" hazerad_file = ",trim(hazerad_file)
+         write(*,*) " hazerad_file = ",trim(hazerad_file)
          write(*,*) "Haze mmr datafile"
          hazemmr_file="None"  ! default file
          call getin("hazemmr_file",hazemmr_file)
-         write(*,*) trim(rname)//" hazemmr_file = ",trim(hazemmr_file)
+         write(*,*) " hazemmr_file = ",trim(hazemmr_file)
          write(*,*) "Haze dens datafile"
          hazedens_file="None"  ! default file
          call getin("hazedens_file",hazedens_file)
-         write(*,*) trim(rname)//" hazedens_file = ",trim(hazedens_file)
+         write(*,*) " hazedens_file = ",trim(hazedens_file)
 
 !***************************************************************
Index: trunk/LMDZ.PLUTO.old/libf/phypluto/optci.F90
===================================================================
--- trunk/LMDZ.PLUTO.old/libf/phypluto/optci.F90	(revision 3357)
+++ trunk/LMDZ.PLUTO.old/libf/phypluto/optci.F90	(revision 3373)
@@ -66,5 +66,6 @@
 
 !     temporary variables for multiple aerosol calculation
-      real*8 atemp, btemp
+      real*8 atemp 
+      real*8 btemp(L_NLAYRAD,L_NSPECTI)
 
 !     variables for k in units m^-1
@@ -186,5 +187,5 @@
 
       DO NW=1,L_NSPECTI
-         DO K=2,L_LEVELS+1
+         DO K=2,L_LEVELS
             do iaer=1,naerkind
                TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
@@ -201,14 +202,15 @@
 
             atemp = 0.
-            btemp = 0.
+            btemp(L,NW) = 0.
+            do iaer=1,naerkind
+               atemp = atemp +                                     &
+                     GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
+                     GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
+               btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
+!     *                    + 1.e-10
+            end do
+            
             if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
-               do iaer=1,naerkind
-                  atemp = atemp +                                     &
-                      GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
-                      GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
-                  btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
-!     *                    + 1.e-10
-               end do
-               WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
+               WBARI(L,nw,ng) = btemp(L,NW)  / DTAUI(L,NW,NG)
             else
                WBARI(L,nw,ng) = 0.0D0
@@ -216,6 +218,6 @@
             endif
 
-            if(btemp .GT. 0.0) then
-               cosbi(L,NW,NG) = atemp/btemp
+            if(btemp(L,NW) .GT. 0.0) then
+               cosbi(L,NW,NG) = atemp/btemp(L,NW)
             else
                cosbi(L,NW,NG) = 0.0D0
@@ -223,29 +225,32 @@
 
          END DO ! L vertical loop
+         ! Last level
          
-     ! Last level
-     
-     L              = L_NLAYRAD
-     K              = 2*L+1
-     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
-
-     atemp = 0.
-     if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
-        do iaer=1,naerkind
-           atemp = atemp + GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
-        end do
-        WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
-     else
-        WBARI(L,nw,ng) = 0.0D0
-        DTAUI(L,NW,NG) = 1.0D-9
-     endif
-
-     if(btemp .GT. 0.0d0) then
-        cosbi(L,NW,NG) = atemp/btemp
-     else
-        cosbi(L,NW,NG) = 0.0D0
-     end if
-
-!     Now the other Gauss points, if needed.
+         L              = L_NLAYRAD
+         K              = 2*L+1
+         DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
+         btemp(L,NW) = 0
+         do iaer=1,naerkind
+               btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,IAER)
+         enddo
+         
+         atemp = 0.
+         if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
+            do iaer=1,naerkind
+               atemp = atemp + GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
+            end do
+            WBARI(L,nw,ng) = btemp(L,NW)  / DTAUI(L,NW,NG)
+         else
+            WBARI(L,nw,ng) = 0.0D0
+            DTAUI(L,NW,NG) = 1.0D-9
+         endif
+
+         if(btemp(L,NW) .GT. 0.0d0) then
+            cosbi(L,NW,NG) = atemp/btemp(L,NW)
+         else
+            cosbi(L,NW,NG) = 0.0D0
+         end if
+
+         !     Now the other Gauss points, if needed.
 
          DO NG=1,L_NGAUSS-1
@@ -256,12 +261,6 @@
                   DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
 
-                  btemp = 0.
                   if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
-
-                     do iaer=1,naerkind
-                        btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
-                     end do
-                     WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
-
+                     WBARI(L,nw,ng) = btemp(L,NW)  / DTAUI(L,NW,NG)
                   else
                      WBARI(L,nw,ng) = 0.0D0
Index: trunk/LMDZ.PLUTO.old/libf/phypluto/physiq.F
===================================================================
--- trunk/LMDZ.PLUTO.old/libf/phypluto/physiq.F	(revision 3357)
+++ trunk/LMDZ.PLUTO.old/libf/phypluto/physiq.F	(revision 3373)
@@ -1163,4 +1163,5 @@
           zdqphot_prec(:,:)=0.
           zdqphot_ch4(:,:)=0.
+          zdqhaze(:,:,:)=0.
           ! Forcing to a fixed haze profile if haze_proffix
           if (haze_proffix.and.i_haze.gt.0.) then