Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 877)
+++ /trunk/LMDZ.GENERIC/README	(revision 878)
@@ -901,4 +901,11 @@
 - Added a keyword to enable ocean runoff in callphys.def (activerunoff)
 
+== 12/02/2013 == JL
+- Follows previous commit by Aymeric about bilinear interpolations:
+   - Extended to all existing continua
+   - generalized bilinearbig to work for various size inputs
+- because N2 and H2O continua databases are smaller, improvement around 15% for
+an earth case.
  
 
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/bilinearbig.F90	(revision 878)
@@ -1,3 +1,3 @@
-      subroutine bilinearbig(x_arr,y_arr,f2d_arr,x_in,y_in,f,a)
+      subroutine bilinearbig(nX,nY,x_arr,y_arr,f2d_arr,x_in,y_in,f,ind)
 
 !     Necessary for interpolation of continuum data
@@ -6,7 +6,5 @@
       implicit none
 
-      integer nX,nY,i,j,a,b
-      parameter(nX=2428)
-      parameter(nY=10)
+      integer nX,nY,i,j,ind,b
 
       real*8 x_in,y_in,x1,x2,y1,y2
@@ -28,9 +26,9 @@
    !! ... and actually calculations only need to be done once
    !! --> Case 1 : we have not calculated yet
-   if ( a == -9999) then
+   if ( ind == -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
+         ind=-1
       else
         i=1
@@ -40,15 +38,15 @@
           i=i+1
           x2=x_arr(i)
-          a=i-1
+          ind=i-1
         end do
       endif
    !! --> case 2 : we already saw we are out of wavenumber range
-   else if ( a == -1) then 
+   else if ( ind == -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)
+      x1=x_arr(ind)
+      x2=x_arr(ind+1)
    endif
 
@@ -74,8 +72,8 @@
       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)
+      f11=f2d_arr(ind,b)
+      f21=f2d_arr(ind+1,b)
+      f12=f2d_arr(ind,b+1)
+      f22=f2d_arr(ind+1,b+1)
 
       call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 878)
@@ -1,3 +1,3 @@
-     subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,a)
+     subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,ind)
 
 !==================================================================
@@ -54,5 +54,5 @@
       integer nres
 
-      integer a
+      integer ind
 
       if(temp.gt.400)then
@@ -109,5 +109,5 @@
       endif
 
-         call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
+         call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
 
          !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90	(revision 878)
@@ -1,3 +1,3 @@
-     subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,a)
+     subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,ind)
 
 !==================================================================
@@ -56,5 +56,5 @@
       integer nres
 
-      integer a
+      integer ind
  
       if(temp.gt.400)then
@@ -112,5 +112,5 @@
       endif
 
-         call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
+         call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
 
          !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_CKD.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_CKD.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_CKD.F90	(revision 878)
@@ -1,3 +1,3 @@
-     subroutine interpolateH2Ocont_CKD(wn,temp,presS,presF,abcoef,firstcall)
+     subroutine interpolateH2Ocont_CKD(wn,temp,presS,presF,abcoef,firstcall,ind)
 
 !==================================================================
@@ -40,5 +40,5 @@
       double precision data_tmp(nT)
 
-      integer k
+      integer k,ind
       logical firstcall
 
@@ -129,131 +129,39 @@
          print*,'   H2O pressure ',presS,' Pa'
          print*,'   air pressure ',presF,' Pa'
+	 
+      endif
 
-         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
-         print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'
+      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS,ind)
+!      print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'
 
-         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
-         print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'
+      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF,ind)
+!      print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'
 
-         print*,'We have ',amagatS,' amagats of H2O vapour'
-         print*,'and ',amagatF,' amagats of air'
+!      print*,'We have ',amagatS,' amagats of H2O vapour'
+!      print*,'and ',amagatF,' amagats of air'
 
-         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
-         abcoef = abcoef*(presS/(presF+presS))      ! take H2O mixing ratio into account
+      abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
+      abcoef = abcoef*(presS/(presF+presS))      ! take H2O mixing ratio into account
                                                     ! abs coeffs are given per molecule of H2O
 
-         Nmolec = (presS+presF)/(kB*temp)           ! assume ideal gas
-         print*,'Total number of molecules per m^3 is',Nmolec
+      Nmolec = (presS+presF)/(kB*temp)           ! assume ideal gas
+!      print*,'Total number of molecules per m^3 is',Nmolec
 
-         abcoef = abcoef*Nmolec/(100.0**2)          ! convert to m^-1
-         print*,'So the total absorption is ',abcoef,' m^-1'
-         print*,'And optical depth / km : ',1000.0*abcoef
+      abcoef = abcoef*Nmolec/(100.0**2)          ! convert to m^-1
+!      print*,'So the total absorption is ',abcoef,' m^-1'
+!      print*,'And optical depth / km : ',1000.0*abcoef
 
+
+      if(wn.gt.500 .and. wn.lt.1400)then
+      elseif(wn.gt.2100 .and. wn.lt.3000)then
       else
+         abcoef = 0.0
+      endif
 
-         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
-         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
+      ! 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.
 
-         abcoef = abcoefS*amagatS + abcoefF*amagatF
-         abcoef = abcoef*(presS/(presF+presS)) 
-
-         Nmolec = (presS+presF)/(kB*temp)
-         abcoef = abcoef*Nmolec/(100.0**2) 
-
-         if(wn.gt.500 .and. wn.lt.1400)then
-         elseif(wn.gt.2100 .and. wn.lt.3000)then
-         else
-            abcoef = 0.0
-         endif
-
-         ! 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 interpolateH2Ocont_CKD
 
-
-!-------------------------------------------------------------------------
-      subroutine bilinearH2Ocont(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=1001)
-      parameter(nY=11)
-
-      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 bilinearH2Ocont:'
-         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)
-      
-!     1st in x-direction
-      fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
-      fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
-      
-!     then in y-direction
-      f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
-
-      return
-    end subroutine bilinearH2Ocont
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateN2N2.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateN2N2.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateN2N2.F90	(revision 878)
@@ -1,3 +1,3 @@
-subroutine interpolateN2N2(wn,temp,pres,abcoef,firstcall)
+subroutine interpolateN2N2(wn,temp,pres,abcoef,firstcall,ind)
 
   !==================================================================
@@ -21,4 +21,5 @@
   double precision temp               ! temperature            (Kelvin)
   double precision pres               ! pressure               (Pascals)
+  integer :: ind
 
   ! output
@@ -104,102 +105,21 @@
      print*,'   pressure ',pres,' Pa'
 
-     call bilinearN2N2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
+  endif
+     call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
 
-     print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
-     print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
+!     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
+!     abcoef=0.
 
-     print*,'We have ',amagat,' amagats of N2'
-     print*,'So the absorption is ',abcoef,' m^-1'
+!     print*,'We have ',amagat,' amagats of N2'
+!     print*,'So the absorption is ',abcoef,' m^-1'
 
-  else
+!     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.
 
-     call bilinearN2N2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
-     abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to 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 interpolateN2N2
 
-
-!-------------------------------------------------------------------------
-      subroutine bilinearN2N2(x_arr,y_arr,f2d_arr,x_in,y_in,f)
-
-      implicit none
-
-      integer nX,nY,i,j,a,b
-      parameter(nX=582)
-      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 bilinearN2N2:'
-         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 bilinearN2N2
Index: /trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 878)
@@ -92,5 +92,5 @@
   !! AS: to save time in computing continuum (see bilinearbig)
   IF (.not.ALLOCATED(indi)) THEN
-      ALLOCATE(indi(L_NSPECTI,ngasmx))
+      ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
       indi = -9999  ! this initial value means "to be calculated"
   ENDIF
@@ -162,12 +162,14 @@
               if(igas.eq.igas_N2)then
 
-                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                 interm = indi(nw,igas,igas)
+                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
+                 indi(nw,igas,igas) = interm
 
               elseif(igas.eq.igas_H2)then
 
                  ! first do self-induced absorption
-                 interm = indi(nw,igas)
+                 interm = indi(nw,igas,igas)
                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
-                 indi(nw,igas) = interm
+                 indi(nw,igas,igas) = interm
 
                  ! then cross-interactions with other gases
@@ -176,9 +178,11 @@
                     dtempc  = 0.0
                     if(jgas.eq.igas_N2)then 
-                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
+                       interm = indi(nw,igas,jgas)
+                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indi(nw,igas,jgas) = interm
                     elseif(jgas.eq.igas_He)then 
-                       interm = indi(nw,jgas)
+                       interm = indi(nw,igas,jgas)
                        call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
-                       indi(nw,jgas) = interm
+                       indi(nw,igas,jgas) = interm
                     endif
                     dtemp = dtemp + dtempc
@@ -191,5 +195,7 @@
                     call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
                  else
-                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+                    interm = indi(nw,igas,igas)
+                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
+                    indi(nw,igas,igas) = interm
                  endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 878)
@@ -87,5 +87,5 @@
   !! AS: to save time in computing continuum (see bilinearbig)
   IF (.not.ALLOCATED(indv)) THEN
-      ALLOCATE(indv(L_NSPECTV,ngasmx))
+      ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
       indv = -9999 ! this initial value means "to be calculated"
   ENDIF
@@ -157,5 +157,7 @@
               if(igas.eq.igas_N2)then
 
-                 !call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                 interm = indv(nw,igas,igas)
+!                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
+                 indv(nw,igas,igas) = interm
                  ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
 
@@ -163,7 +165,7 @@
 
                  ! first do self-induced absorption
-                 interm = indv(nw,igas)
+                 interm = indv(nw,igas,igas)
                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
-                 indv(nw,igas) = interm
+                 indv(nw,igas,igas) = interm
 
                  ! then cross-interactions with other gases
@@ -172,10 +174,12 @@
                     dtempc  = 0.0
                     if(jgas.eq.igas_N2)then 
-                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
+                       interm = indv(nw,igas,jgas)
+                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indv(nw,igas,jgas) = interm
                        ! should be irrelevant in the visible
                     elseif(jgas.eq.igas_He)then 
-                       interm = indv(nw,jgas)
+                       interm = indv(nw,igas,jgas)
                        call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
-                       indv(nw,jgas) = interm
+                       indv(nw,igas,jgas) = interm
                     endif
                     dtemp = dtemp + dtempc
@@ -188,5 +192,7 @@
                     call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
                  else
-                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+                    interm = indv(nw,igas,igas)
+                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
+                    indv(nw,igas,igas) = interm
                  endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 877)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 878)
@@ -67,6 +67,6 @@
 
       !! AS: introduced to avoid doing same computations again for continuum
-      INTEGER, DIMENSION(:,:), ALLOCATABLE :: indi
-      INTEGER, DIMENSION(:,:), ALLOCATABLE :: indv
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indi
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indv
 
       !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011  
