Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 872)
+++ /trunk/LMDZ.GENERIC/README	(revision 873)
@@ -878,2 +878,20 @@
 - Users can still use e.g. H2_ but H2 also works
 - Simplified condense_cloud so that igas_CO2 is used directly
+
+== 09/02/2013 == AS
+- Optimized calculations for continuum (done for H2 and He, to be done for others)
+      - new common bilinear interpolation routine (bilinearbig)
+      - optimization: only one calculation is actually needed
+                      to find indexes of wavelength for bilinear interpolation
+                      ... because this will not change with level and integration step!
+      - optimization: use while loop in bilinearbig
+      - completely similar results obtained (test case for a gas giant, many simulated days)
+  NB: those changes really improve gcm speed (factor 2.2 for whole model!)
+      continuum was very expensive, now very cheap
+      --> e.g. 1 day, 25 dyn ts, 5 phys ts
+      --> before: 243 seconds (including 120 seconds for continuum bilinear interpolation)
+      --> after: 108 seconds
+- Corrected a bug: Continuum in inifis instead of continuum
+    ... until now, most users (unbeknownst to them) were running with the continuum by default!
+- Cosmetic changes in optcv (mostly spaces and line breaks) 
+    ... so that comparisons with optci are easy e.g. through vimdiff
Index: /trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90	(revision 873)
+++ /trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90	(revision 873)
@@ -0,0 +1,84 @@
+      subroutine bilinearbig(x_arr,y_arr,f2d_arr,x_in,y_in,f,a)
+
+!     Necessary for interpolation of continuum data
+!     optimized by A. Spiga 01/2013 
+
+      implicit none
+
+      integer nX,nY,i,j,a,b
+      parameter(nX=2428)
+      parameter(nY=10)
+
+      real*8 x_in,y_in,x1,x2,y1,y2
+      real*8 f,f11,f12,f21,f22,fA,fB
+      real*8 x_arr(nX)
+      real*8 y_arr(nY)
+      real*8 f2d_arr(nX,nY)
+      real*8,save :: x,y
+
+      integer strlen
+      character*100 label
+      label='subroutine bilinear'
+
+      x=x_in
+      y=y_in
+
+
+   !! AS: important to optimize here because the array is quite large
+   !! ... and actually calculations only need to be done once
+   !! --> Case 1 : we have not calculated yet
+   if ( a == -9999) then
+      !1st check we're within the wavenumber range
+      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
+         f=0.0D+0
+         a=-1
+      else
+        i=1
+        x2=x_arr(i)
+        do while ( x2 .le. x )
+          x1=x2
+          i=i+1
+          x2=x_arr(i)
+          a=i-1
+        end do
+      endif
+   !! --> case 2 : we already saw we are out of wavenumber range
+   else if ( a == -1) then 
+      f=0.0D+0
+      return
+   !! --> case 3 : we already determined a -- so we just proceed
+   else
+      x1=x_arr(a)
+      x2=x_arr(a+1)
+   endif
+
+!     ... and for y within the temperature range
+      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
+         write(*,*) 'Warning from bilinearH2H2:'
+         write(*,*) 'Outside continuum temperature range!'
+         if(y.lt.y_arr(1))then
+            y=y_arr(1)+0.01
+         endif
+         if(y.gt.y_arr(nY))then
+            y=y_arr(nY)-0.01
+         endif
+      else
+        j=1
+        y2=y_arr(j)
+        do while ( y2 .le. y )
+          y1=y2
+          j=j+1
+          y2=y_arr(j)
+          b=j-1
+        end do
+      endif
+
+      f11=f2d_arr(a,b)
+      f21=f2d_arr(a+1,b)
+      f12=f2d_arr(a,b+1)
+      f22=f2d_arr(a+1,b+1)
+
+      call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
+
+      return
+    end subroutine bilinearbig
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 873)
@@ -8,5 +8,5 @@
      &   , iaervar,iddist,topdustref,callstats,calleofdump              &
      &   , enertest                                                     &
-     &   , callgasvis,Continuum,H2Ocont_simple,graybody                 &
+     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
      &   , Nmix_co2, radfixed, dusttau                                  &
      &   , meanOLR, specOLR                                             &
@@ -35,5 +35,5 @@
      &   , season,diurnal,tlocked,lwrite                                &
      &   , callstats,calleofdump                                        &
-     &   , callgasvis,Continuum,H2Ocont_simple,graybody
+     &   , callgasvis,continuum,H2Ocont_simple,graybody
 
       logical enertest
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 873)
@@ -203,7 +203,7 @@
          write(*,*) "call continuum opacities in radiative transfer ?",
      &              "(matters only if callrad=T)"
-         Continuum=.true. ! default value
-         call getin("Continuum",Continuum)
-         write(*,*) " Continuum = ",Continuum
+         continuum=.true. ! default value
+         call getin("continuum",continuum)
+         write(*,*) " continuum = ",continuum
 
          write(*,*) "use analytic function for H2O continuum ?"
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 873)
@@ -1,3 +1,3 @@
-     subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall)
+     subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,a)
 
 !==================================================================
@@ -15,4 +15,5 @@
 
       use datafile_mod, only: datadir
+
       implicit none
 
@@ -52,4 +53,6 @@
       double precision blah, Ttemp
       integer nres
+
+      integer a
 
       if(temp.gt.400)then
@@ -104,124 +107,19 @@
          print*,'   pressure ',pres,' Pa'
 
-         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
+      endif
 
-         print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
-         print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
+         call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
+
+         !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
+         !print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
 
          abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
 
-         print*,'We have ',amagat,' amagats of H2'
-         print*,'So the absorption is ',abcoef,' m^-1'
-
-               !open(117,file='T_array.dat')
-               !do iT=1,nT
-               !   write(117,*), temp_arr(iT)
-               !end do
-               !close(117)
-
-               !open(115,file='nu_array.dat')
-               !do k=1,nS
-               !   write(115,*), wn_arr(k)
-               !end do
-               !close(115)
-
-               !open(113,file='abs_array.dat')
-               !do iT=1,nT
-               !   do k=1,nS
-               !      write(113,*), abs_arr(k,iT)*losch**2
-               !   end do
-               !end do
-               !close(113)
-
-      else
-  
-         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
-         abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
+         !print*,'We have ',amagat,' amagats of H2'
+         !print*,'So the absorption is ',abcoef,' m^-1'
 
          ! unlike for Rayleigh scattering, we do not currently weight by the BB function
          ! however our bands are normally thin, so this is no big deal.
 
-      endif
-
       return
     end subroutine interpolateH2H2
-
-
-!-------------------------------------------------------------------------
-      subroutine bilinearH2H2(x_arr,y_arr,f2d_arr,x_in,y_in,f)
-!     Necessary for interpolation of continuum data
-
-      implicit none
-
-      integer nX,nY,i,j,a,b
-      parameter(nX=2428)
-      parameter(nY=10)
-
-      real*8 x_in,y_in,x,y,x1,x2,y1,y2
-      real*8 f,f11,f12,f21,f22,fA,fB
-      real*8 x_arr(nX)
-      real*8 y_arr(nY)
-      real*8 f2d_arr(nX,nY)
-      
-      integer strlen
-      character*100 label
-      label='subroutine bilinear'
-
-      x=x_in
-      y=y_in
-
-
-!     1st check we're within the wavenumber range
-      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
-         f=0.0D+0
-         return
-      else
-         
-!     in the x (wavenumber) direction 1st
-         i=1
- 10      if (i.lt.(nX+1)) then
-            if (x_arr(i).gt.x) then
-               x1=x_arr(i-1)
-               x2=x_arr(i)
-               a=i-1
-               i=9999
-            endif
-            i=i+1
-            goto 10
-         endif
-      endif
-      
-      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
-         write(*,*) 'Warning from bilinearH2H2:'
-         write(*,*) 'Outside continuum temperature range!'
-         if(y.lt.y_arr(1))then
-            y=y_arr(1)+0.01
-         endif
-         if(y.gt.y_arr(nY))then
-            y=y_arr(nY)-0.01
-         endif
-      else
-
-!     in the y (temperature) direction 2nd
-         j=1
- 20      if (j.lt.(nY+1)) then
-            if (y_arr(j).gt.y) then
-               y1=y_arr(j-1)
-               y2=y_arr(j)
-               b=j-1
-               j=9999
-            endif
-            j=j+1
-            goto 20
-         endif
-      endif
-      
-      f11=f2d_arr(a,b)
-      f21=f2d_arr(a+1,b)
-      f12=f2d_arr(a,b+1)
-      f22=f2d_arr(a+1,b+1)
-
-      call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
-
-      return
-    end subroutine bilinearH2H2
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90	(revision 873)
@@ -1,3 +1,3 @@
-     subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall)
+     subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,a)
 
 !==================================================================
@@ -15,4 +15,5 @@
 
       use datafile_mod, only: datadir
+
       implicit none
 
@@ -55,5 +56,6 @@
       integer nres
 
-
+      integer a
+ 
       if(temp.gt.400)then
          print*,'Your temperatures are too high for this H2-He CIA dataset.'
@@ -108,105 +110,20 @@
          print*,'   and He partial pressure     ',presHe,' Pa'
 
-         call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
+      endif
 
-         print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
-         print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
+         call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
+
+         !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
+         !print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
 
          abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1
 
-         print*,'We have ',amagatH2,' amagats of H2'
-         print*,'and     ',amagatHe,' amagats of He'
-         print*,'So the absorption is ',abcoef,' m^-1'
-
-      else
-  
-         call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
-         abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1
+         !print*,'We have ',amagatH2,' amagats of H2'
+         !print*,'and     ',amagatHe,' amagats of He'
+         !print*,'So the absorption is ',abcoef,' m^-1'
 
          ! unlike for Rayleigh scattering, we do not currently weight by the BB function
          ! however our bands are normally thin, so this is no big deal.
 
-      endif
-
       return
     end subroutine interpolateH2He
-
-
-!-------------------------------------------------------------------------
-      subroutine bilinearH2He(x_arr,y_arr,f2d_arr,x_in,y_in,f)
-!     Necessary for interpolation of continuum data
-
-      implicit none
-
-      integer nX,nY,i,j,a,b
-      parameter(nX=2428)
-      parameter(nY=10)
-
-      real*8 x_in,y_in,x,y,x1,x2,y1,y2
-      real*8 f,f11,f12,f21,f22,fA,fB
-      real*8 x_arr(nX)
-      real*8 y_arr(nY)
-      real*8 f2d_arr(nX,nY)
-      
-      integer strlen
-      character*100 label
-      label='subroutine bilinear'
-
-      x=x_in
-      y=y_in
-
-!     1st check we're within the wavenumber range
-      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
-         f=0.0D+0
-         return
-      else
-         
-!     in the x (wavenumber) direction 1st
-         i=1
- 10      if (i.lt.(nX+1)) then
-            if (x_arr(i).gt.x) then
-               x1=x_arr(i-1)
-               x2=x_arr(i)
-               a=i-1
-               i=9999
-            endif
-            i=i+1
-            goto 10
-         endif
-      endif
-      
-      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
-         write(*,*) 'Warning from bilinearH2He:'
-         write(*,*) 'Outside continuum temperature range!'
-         if(y.lt.y_arr(1))then
-            y=y_arr(1)+0.01
-         endif
-         if(y.gt.y_arr(nY))then
-            y=y_arr(nY)-0.01
-         endif
-      else
-
-!     in the y (temperature) direction 2nd
-         j=1
- 20      if (j.lt.(nY+1)) then
-            if (y_arr(j).gt.y) then
-               y1=y_arr(j-1)
-               y2=y_arr(j)
-               b=j-1
-               j=9999
-            endif
-            j=j+1
-            goto 20
-         endif
-      endif
-      
-      f11=f2d_arr(a,b)
-      f21=f2d_arr(a+1,b)
-      f12=f2d_arr(a,b+1)
-      f22=f2d_arr(a+1,b+1)
-
-      call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
-
-      
-      return
-    end subroutine bilinearH2He
Index: /trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 873)
@@ -4,5 +4,5 @@
 
   use radinc_h
-  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep
+  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi
   use gases_h
   implicit none
@@ -70,7 +70,10 @@
 
   ! variables for k in units m^-1
-  real*8 rho, dz(L_LEVELS)
+  real*8 dz(L_LEVELS)
+  !real*8 rho !! see test below
 
   integer igas, jgas
+
+  integer interm
 
   !--- Kasting's CIA ----------------------------------------
@@ -87,4 +90,10 @@
   !----------------------------------------------------------
 
+  !! AS: to save time in computing continuum (see bilinearbig)
+  IF (.not.ALLOCATED(indi)) THEN
+      ALLOCATE(indi(L_NSPECTI,ngasmx))
+      indi = -9999  ! this initial value means "to be calculated"
+  ENDIF
+
   !=======================================================================
   !     Determine the total gas opacity throughout the column, for each
@@ -133,9 +142,10 @@
 
   do K=2,L_LEVELS
-     do nw=1,L_NSPECTI
+
+     do NW=1,L_NSPECTI
 
         DCONT = 0.0 ! continuum absorption
 
-        if(Continuum.and.(.not.graybody))then
+        if(continuum.and.(.not.graybody))then
            ! include continua if necessary
            wn_cont = dble(wnoi(nw))
@@ -157,5 +167,7 @@
 
                  ! first do self-induced absorption
-                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                 interm = indi(nw,igas)
+                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
+                 indi(nw,igas) = interm
 
                  ! then cross-interactions with other gases
@@ -166,5 +178,7 @@
                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
                     elseif(jgas.eq.igas_He)then 
-                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
+                       interm = indi(nw,jgas)
+                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indi(nw,jgas) = interm
                     endif
                     dtemp = dtemp + dtempc
@@ -228,5 +242,5 @@
 
               KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*     &
-                   (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                   & 
+                   (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                   &
                    GASI(MT(K),MP(K),NVAR(K),NW,NG))
 
@@ -242,4 +256,5 @@
                    (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                 &
                    GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
+
            endif
 
@@ -252,5 +267,6 @@
 
            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
-           DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
+           DTAUKI(K,nw,ng) = TAUGAS & 
+                             + DCONT ! For parameterized continuum absorption
 
            do iaer=1,naerkind
@@ -278,11 +294,11 @@
   !     we need to calculate the scattering albedo and asymmetry factors
 
+  do iaer=1,naerkind
   DO NW=1,L_NSPECTI
      DO K=2,L_LEVELS+1
-        do iaer=1,naerkind
            TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
-        end do
      ENDDO
   ENDDO
+  end do
 
   DO NW=1,L_NSPECTI
Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 873)
@@ -4,5 +4,5 @@
 
   use radinc_h
-  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
+  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv
   use gases_h
 
@@ -35,6 +35,7 @@
   !-------------------------------------------------------------------
 
+#include "comcstfi.h"
 #include "callkeys.h"
-#include "comcstfi.h"
+
 
   real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
@@ -46,26 +47,30 @@
   real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
   real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-  real*8 TAURAY(L_NSPECTV)
 
   ! for aerosols
-  real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
-  real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
-  real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
-  real*8 TAUAERO(L_LEVELS+1,NAERKIND)
-  real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
-  real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
-
-  integer L, NW, NG, K, NG1(L_NSPECTV), LK, IAER
+  real*8  QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
+  real*8  QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
+  real*8  GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
+  real*8  TAUAERO(L_LEVELS+1,NAERKIND)
+  real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
+  real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
+
+  integer L, NW, NG, K, LK, IAER
   integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
   real*8  ANS, TAUGAS
+  real*8  TAURAY(L_NSPECTV)
   real*8  TRAY(L_LEVELS,L_NSPECTV)
+  real*8  TRAYAER
   real*8  DPR(L_LEVELS), U(L_LEVELS)
   real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
 
-  real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
+  real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
+  real*8 DCONT
+  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
+  double precision p_cross
 
   ! variable species mixing ratio variables
-  real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
-  real*8 KCOEF(4)
+  real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
+  real*8  KCOEF(4)
   integer NVAR(L_LEVELS)
 
@@ -74,9 +79,15 @@
 
   ! variables for k in units m^-1
-  double precision wn_cont, p_cont, p_air, T_cont, dtemp
-  double precision p_cross
-  real*8 dz(L_LEVELS), DCONT
+  real*8 dz(L_LEVELS)
 
   integer igas, jgas
+
+  integer interm
+
+  !! AS: to save time in computing continuum (see bilinearbig)
+  IF (.not.ALLOCATED(indv)) THEN
+      ALLOCATE(indv(L_NSPECTV,ngasmx))
+      indv = -9999 ! this initial value means "to be calculated"
+  ENDIF
 
   !=======================================================================
@@ -95,5 +106,5 @@
      ! if we have continuum opacities, we need dz
      if(kastprof)then
-        dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
+        dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K))
         U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
      else
@@ -109,21 +120,20 @@
      end do
 
+
      DO NW=1,L_NSPECTV
+        do iaer=1,naerkind
+           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
+        end do
         TRAY(K,NW)   = TAURAY(NW) * DPR(K)
-
-        do iaer=1,naerkind
-           TAEROS(K,NW,IAER) = TAUAERO(K,IAER)  * QXVAER(K,NW,IAER)
-        end do
-
      END DO
-  end do
-
-  !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
+  end do                    ! levels
 
   !     we ignore K=1...
   do K=2,L_LEVELS
+
      do NW=1,L_NSPECTV
 
         TRAYAER = TRAY(K,NW)
+        !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
         do iaer=1,naerkind
            TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
@@ -132,5 +142,5 @@
         DCONT = 0.0 ! continuum absorption
 
-        if(callgasvis.and.continuum.and.(.not.graybody))then
+        if(continuum.and.(.not.graybody).and.callgasvis)then
            ! include continua if necessary
            wn_cont = dble(wnov(nw))
@@ -153,15 +163,21 @@
 
                  ! first do self-induced absorption
-                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                 interm = indv(nw,igas)
+                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
+                 indv(nw,igas) = interm
 
                  ! then cross-interactions with other gases
                  do jgas=1,ngasmx
                     p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
-                    if(jgas.eq.igas_N2)then
-                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
+                    dtempc  = 0.0
+                    if(jgas.eq.igas_N2)then 
+                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
                        ! should be irrelevant in the visible
                     elseif(jgas.eq.igas_He)then 
-                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
+                       interm = indv(nw,jgas)
+                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indv(nw,jgas) = interm
                     endif
+                    dtemp = dtemp + dtempc
                  enddo
 
@@ -181,16 +197,17 @@
            enddo
 
-           DCONT = DCONT*dz(k) 
+           DCONT = DCONT*dz(k)
+
         endif
 
-        do NG=1,L_NGAUSS-1
-
-           !=======================================================================
-           !     Now compute TAUGAS
-           !     Interpolate between water mixing ratios
-           !     WRATIO = 0.0 if the requested water amount is equal to, or outside the
-           !     the water data range
-
-           if (L_REFVAR.eq.1)then ! added by RW for special no variable case
+        do ng=1,L_NGAUSS-1
+
+           ! Now compute TAUGAS
+
+           ! Interpolate between water mixing ratios
+           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
+           ! the water data range
+
+           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
               KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
               KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
@@ -198,4 +215,5 @@
               KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
            else
+
               KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*    &
                    (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -                  &
@@ -213,29 +231,29 @@
                    (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                &
                    GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
+
            endif
 
-           !     Interpolate the gaseous k-coefficients to the requested T,P values
-
-           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +           &
+           ! Interpolate the gaseous k-coefficients to the requested T,P values
+
+           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
                 LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
 
-           TAUGAS = U(k)*ANS
-
-           !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
+           TAUGAS  = U(k)*ANS
+
            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
-           DTAUKV(K,nw,ng) = TAUGAS + TRAYAER  & ! TRAYAER includes all scattering contributions
-                + DCONT             ! for continuum absorption
+           DTAUKV(K,nw,ng) = TAUGAS & 
+                             + TRAYAER & ! TRAYAER includes all scattering contributions
+                             + DCONT ! For parameterized continuum aborption
 
         end do
 
-
-        !     Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
-        !     which holds continuum opacity only
-
-        NG = L_NGAUSS
+        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
+        ! which holds continuum opacity only
+
+        NG              = L_NGAUSS
         DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
+
         do iaer=1,naerkind
-           DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
-           !     &                           + DCONT          ! For parameterized continuum absorption
+           DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) +  TAEROS(K,NW,IAER)
         end do ! a bug was here!
 
@@ -248,37 +266,33 @@
   !     we need to calculate the scattering albedo and asymmetry factors
 
+  do iaer=1,naerkind
   DO NW=1,L_NSPECTV
-     DO K=2,L_LEVELS
-        do iaer=1,naerkind
+     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)
-        end do
      ENDDO
   ENDDO
-
+  end do
 
   DO NW=1,L_NSPECTV
-     DO NG=1,L_NGAUSS
-        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 NG=1,L_NGAUSS
+     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)
+              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/DTAUV(L,nw,ng)
-
-        END DO
+
+      END DO ! L vertical loop
 
         !     No vertical averaging on bottom layer
@@ -299,12 +313,10 @@
         WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
 
-     END DO                 ! NG gauss point loop
+     END DO                 ! NG Gauss loop
   END DO                    ! NW spectral loop
 
-
-
   ! Total extinction optical depths
 
-  DO NW=1,L_NSPECTV        
+  DO NW=1,L_NSPECTV       
      DO NG=1,L_NGAUSS       ! full gauss loop
         TAUV(1,NW,NG)=0.0D0
@@ -321,4 +333,6 @@
 
 
-  RETURN
-END SUBROUTINE OPTCV
+  return
+
+
+end subroutine optcv
Index: /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 873)
@@ -66,4 +66,8 @@
       REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90
 
+      !! AS: introduced to avoid doing same computations again for continuum
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: indi
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: indv
+
       !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011  
       REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv
Index: /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 872)
+++ /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 873)
@@ -55,4 +55,6 @@
 
       double precision testcont ! for continuum absorption initialisation
+
+      integer :: dummy
 
 !=======================================================================
@@ -596,5 +598,6 @@
 
             ! first do self-induced absorption
-            call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.)
+            dummy = -9999
+            call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
             ! then cross-interactions with other gases
             do jgas=1,ngasmx
@@ -602,5 +605,6 @@
                   call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.)
                elseif (jgas .eq. igas_He) then
-                  call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.)
+                  dummy = -9999
+                  call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
                endif
             enddo
Index: /trunk/LMDZ.UNIVERSAL/README.aymeric
===================================================================
--- /trunk/LMDZ.UNIVERSAL/README.aymeric	(revision 872)
+++ /trunk/LMDZ.UNIVERSAL/README.aymeric	(revision 873)
@@ -32,4 +32,10 @@
 makelmdz_fcm -d 256x256x20 -arch GNOMEp_large -parallel mpi -p generic gcm
 --------------------------------------------------
+
+profiling
+ajouter -pg
+compiler lancer
+puis gprof executable > prof.out
+
 
 --------------------------------------------------
Index: /trunk/LMDZ.UNIVERSAL/libf/phygeneric/bilinearbig.F90
===================================================================
--- /trunk/LMDZ.UNIVERSAL/libf/phygeneric/bilinearbig.F90	(revision 873)
+++ /trunk/LMDZ.UNIVERSAL/libf/phygeneric/bilinearbig.F90	(revision 873)
@@ -0,0 +1,1 @@
+link ../../../LMDZ.GENERIC/libf/phystd/bilinearbig.F90
