Index: trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90	(revision 907)
+++ trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90	(revision 918)
@@ -19,4 +19,5 @@
       label='subroutine bilinear'
 
+
       x=x_in
       y=y_in
@@ -29,4 +30,5 @@
       if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
          ind=-1
+         return
       else
         i=1
@@ -72,5 +74,5 @@
         end do
       endif
-
+      
       f11=f2d_arr(ind,b)
       f21=f2d_arr(ind+1,b)
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 907)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 918)
@@ -223,7 +223,7 @@
          banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
 
-         call sugas_corrk       ! set up gaseous absorption properties
          call setspi            ! basic infrared properties
          call setspv            ! basic visible properties
+         call sugas_corrk       ! set up gaseous absorption properties
          call suaer_corrk       ! set up aerosol optical properties
 
Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 907)
+++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 918)
@@ -57,5 +57,5 @@
 
   real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
-  real*8 DCONT
+  real*8 DCONT,DAERO
   double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
   double precision p_cross
@@ -67,5 +67,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
@@ -132,18 +133,21 @@
         LKCOEF(K,LK) = LCOEF(LK)
      end do
-
-
+  end do                    ! levels
+
+
+  do iaer=1,naerkind
      DO NW=1,L_NSPECTI
-        do iaer=1,naerkind
+        do K=2,L_LEVELS
            TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
-        end do
+        end do                    ! levels
      END DO
-  end do                    ! levels
-
-  do K=2,L_LEVELS
-
-     do NW=1,L_NSPECTI
-
-        DCONT = 0.0 ! continuum absorption
+  end do
+
+  do NW=1,L_NSPECTI
+
+     do K=2,L_LEVELS
+
+! continuum absorption
+        DCONT = 0.0 
 
         if(continuum.and.(.not.graybody))then
@@ -221,14 +225,6 @@
         endif
 
-        ! RW 7/3/12: already done above
-        !if(.not.Continuum)then
-        !   DCONT=0.0
-        !endif
-
-        !--- Kasting's CIA ----------------------------------------
-        !DCO2   = dz(k)*Ci(nw)*(1.2859*PMID(k)/1000.0)*(TMID(k)/300.)**Ti(nw)
-        !DCO2 = 130*Ci(nw)*(pmid(k)/1013.25)**2*(tmid(k)/300.)**Ti(nw) * dz(k)
-        ! these two have been verified to give the same results
-        !----------------------------------------------------------
+! aerosol absorption
+	DAERO=SUM(TAEROS(K,NW,1:naerkind))
 
         do ng=1,L_NGAUSS-1
@@ -273,10 +269,7 @@
 
            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
-           DTAUKI(K,nw,ng) = TAUGAS & 
-                             + DCONT ! For parameterized continuum absorption
-
-           do iaer=1,naerkind
-              DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
-           end do ! a bug was here!
+           DTAUKI(K,nw,ng) = TAUGAS    & 
+                             + DCONT   & ! For parameterized continuum absorption
+			     + DAERO     ! For aerosol absorption
 
         end do
@@ -286,9 +279,7 @@
 
         NG              = L_NGAUSS
-        DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
-
-        do iaer=1,naerkind
-           DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) +  TAEROS(K,NW,IAER)
-        end do ! a bug was here!
+        DTAUKI(K,nw,ng) = 0.d0      & 
+                          + DCONT   & ! For parameterized continuum absorption
+	                  + DAERO     ! For aerosol absorption
 
      end do
@@ -301,10 +292,18 @@
 
   do iaer=1,naerkind
-  DO NW=1,L_NSPECTI
+    DO NW=1,L_NSPECTI
      DO K=2,L_LEVELS+1
            TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
      ENDDO
-  ENDDO
+    ENDDO
   end do
+  
+  DO NW=1,L_NSPECTI
+     DO L=1,L_NLAYRAD
+        K              = 2*L+1
+        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
+     END DO ! L vertical loop
+  END DO                    ! NW spectral loop
+  
 
   DO NW=1,L_NSPECTI
@@ -316,5 +315,4 @@
 
         atemp = 0.
-        btemp = 0.
         if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
            do iaer=1,naerkind
@@ -322,8 +320,6 @@
                    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
@@ -331,6 +327,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
@@ -348,11 +344,7 @@
               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
@@ -370,11 +362,6 @@
   ! Total extinction optical depths
 
-  DO NW=1,L_NSPECTI       
-     DO NG=1,L_NGAUSS       ! full gauss loop
-        TAUI(1,NW,NG)=0.0D0
-        DO L=1,L_NLAYRAD
-           TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
-        END DO
-
+  DO NG=1,L_NGAUSS       ! full gauss loop
+     DO NW=1,L_NSPECTI       
         TAUCUMI(1,NW,NG)=0.0D0
         DO K=2,L_LEVELS
@@ -395,4 +382,11 @@
   !enddo
   !close(127)
+  
+!  print*,'WBARI'
+!  print*,WBARI
+!  print*,'DTAUI'
+!  print*,DTAUI
+!  call abort
+  
 
   return
Index: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 907)
+++ trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 918)
@@ -66,5 +66,5 @@
 
   real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
-  real*8 DCONT
+  real*8 DCONT,DAERO
   double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
   double precision p_cross
@@ -76,5 +76,7 @@
 
   ! temporary variables for multiple aerosol calculation
-  real*8 atemp, btemp, ctemp
+  real*8 atemp(L_NLAYRAD,L_NSPECTV)
+  real*8 btemp(L_NLAYRAD,L_NSPECTV)
+  real*8 ctemp(L_NLAYRAD,L_NSPECTV)
 
   ! variables for k in units m^-1
@@ -119,14 +121,20 @@
         LKCOEF(K,LK) = LCOEF(LK)
      end do
-
-
-     DO NW=1,L_NSPECTV
-        do iaer=1,naerkind
+  end do                    ! levels
+
+
+  do iaer=1,naerkind
+     do NW=1,L_NSPECTV
+        do K=2,L_LEVELS
            TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
-        end do
+        end do                    ! levels
+     end do
+  end do
+  do NW=1,L_NSPECTV
+     do K=2,L_LEVELS
         TRAY(K,NW)   = TAURAY(NW) * DPR(K)
-     END DO
-  end do                    ! levels
-
+     end do                    ! levels
+  end do
+  
   !     we ignore K=1...
   do K=2,L_LEVELS
@@ -273,30 +281,29 @@
 
   do iaer=1,naerkind
+    DO NW=1,L_NSPECTV
+      DO K=2,L_LEVELS   ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
+           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
+      ENDDO
+    ENDDO
+  end do
+
   DO NW=1,L_NSPECTV
-     DO K=2,L_LEVELS   ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
-           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
-     ENDDO
-  ENDDO
-  end do
-
-  DO NW=1,L_NSPECTV
-    DO NG=1,L_NGAUSS
+     DO L=1,L_NLAYRAD
+        K              = 2*L+1
+	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind))
+        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
+	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
+	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
+	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
+     END DO ! L vertical loop
+  END DO                    ! NW spectral loop
+
+  DO NG=1,L_NGAUSS
+    DO NW=1,L_NSPECTV
      DO L=1,L_NLAYRAD-1
 
         K              = 2*L+1
         DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
-
-        atemp = 0.
-        btemp = TRAY(K,NW) + TRAY(K+1,NW)
-        ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
-           do iaer=1,naerkind
-              atemp = atemp +                                     &
-                   GVAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
-                   GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
-              btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
-              ctemp = ctemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
-           end do
-           WBARV(L,nw,ng) = ctemp / DTAUV(L,nw,ng)
-           COSBV(L,NW,NG) = atemp/btemp
+        WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
 
       END DO ! L vertical loop
@@ -307,23 +314,12 @@
         K              = 2*L+1
         DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
-
-        atemp=0.
-        btemp=TRAY(K,NW)
-        ctemp=0.9999*TRAY(K,NW)
-        do iaer=1,naerkind
-           atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) 
-           btemp = btemp + TAUAEROLK(K,NW,IAER) 
-           ctemp = ctemp + TAUAEROLK(K,NW,IAER)
-        end do
-        COSBV(L,NW,NG) = atemp/btemp
-        WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
-
-     END DO                 ! NG Gauss loop
-  END DO                    ! NW spectral loop
+        WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
+     END DO                 ! NW spectral loop
+  END DO                    ! NG Gauss loop
 
   ! Total extinction optical depths
 
-  DO NW=1,L_NSPECTV       
-     DO NG=1,L_NGAUSS       ! full gauss loop
+  DO NG=1,L_NGAUSS       ! full gauss loop
+     DO NW=1,L_NSPECTV       
         TAUV(1,NW,NG)=0.0D0
         DO L=1,L_NLAYRAD
@@ -335,6 +331,6 @@
            TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
         END DO
-     END DO                 ! end full gauss loop
-  END DO
+     END DO            
+  END DO                 ! end full gauss loop
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 907)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 918)
@@ -519,5 +519,5 @@
          else
             print*,'WARNING! Thermal conduction in the soil turned off'
-            capcal(:)=1.e16 !1.e6
+            capcal(:)=1.e6
             fluxgrd(:)=0.
             if(noradsurf)then
@@ -1697,6 +1697,6 @@
 	  endif
 	  if (corrk) then
-!	     call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
-!	     call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
+	     call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
+	     call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
 	  endif 
           if(watercond) then
@@ -1741,7 +1741,7 @@
 
           if(watercond.or.CLFvarying)then
-!             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
-!             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
-!             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
+             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
+             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
+             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
              call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
           endif
