Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 716)
@@ -235,7 +235,10 @@
                           (  0.75 * QREFvis3d(ig,l,iaer) /        &
                           ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
-                            pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
-                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the 
-                            CLF1				    !   clear=false/pq=0 case
+                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
+                          !( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the 
+                          !  CLF1				    !   clear=false/pq=0 case
+                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
+                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   &  
+                            CLF1
 
                      aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0)) 
Index: trunk/LMDZ.GENERIC/libf/phystd/bilinear.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/bilinear.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/bilinear.F90	(revision 716)
@@ -1,80 +1,18 @@
 !-------------------------------------------------------------------------
-      subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x_in,y_in,f)
-!     Necessary for interpolation of continuum data
+subroutine bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)  
+! Used for interpolation of continuum data
 
-      implicit none
+  implicit none
 
-      integer nX,nY,i,j,a,b
+  real*8 x,y,x1,x2,y1,y2
+  real*8 f,f11,f12,f21,f22,fA,fB
 
-      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'
+  ! 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)
 
-      x=x_in
-      y=y_in
+  ! then in y-direction
+  f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
 
-!     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 bilinear.for:'
-         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 bilinear
+  return
+end subroutine bilinear
Index: trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 716)
@@ -40,4 +40,5 @@
             ! all values at 300 K from Engineering Toolbox
             if(gnom(igas).eq.'CO2')then
+               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for Mars conditions) 
                cpp_c   = cpp_c   + 0.846*gfrac(igas)
                mugaz_c = mugaz_c + 44.01*gfrac(igas)
@@ -48,7 +49,16 @@
                cpp_c   = cpp_c   + 14.31*gfrac(igas)
                mugaz_c = mugaz_c + 2.01*gfrac(igas)
+            elseif(gnom(igas).eq.'He_')then
+               cpp_c   = cpp_c   + 5.19*gfrac(igas)
+               mugaz_c = mugaz_c + 4.003*gfrac(igas)
             elseif(gnom(igas).eq.'H2O')then
                cpp_c   = cpp_c   + 1.864*gfrac(igas)
                mugaz_c = mugaz_c + 18.02*gfrac(igas)
+            elseif(gnom(igas).eq.'SO2')then
+               cpp_c   = cpp_c   + 0.64*gfrac(igas)
+               mugaz_c = mugaz_c + 64.066*gfrac(igas)
+            elseif(gnom(igas).eq.'H2S')then
+               cpp_c   = cpp_c   + 1.003*gfrac(igas) ! from wikipedia...
+               mugaz_c = mugaz_c + 34.08*gfrac(igas)
             elseif(gnom(igas).eq.'CH4')then
                cpp_c   = cpp_c   + 2.226*gfrac(igas)
Index: trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 716)
@@ -74,4 +74,6 @@
             elseif(gnom(igas).eq.'He_')then
                print*,'Helium not ready yet!'
+               tauconsti(igas) = 0.0
+               call abort
             else
                print*,'Error in calc_rayleigh: Gas species not recognised!'
@@ -117,4 +119,7 @@
                   elseif(gnom(igas).eq.'H2_')then
                      tauvari(igas) = 1.0/wl**4 
+                  elseif(gnom(igas).eq.'He_')then
+                     print*,'Helium not ready yet!'
+                     tauvari(igas) = 0.0
                   else
                      print*,'Error in calc_rayleigh: Gas species not recognised!'
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 716)
@@ -7,5 +7,4 @@
           reffrad,tau_col,cloudfrac,totcloudfrac,              &
           clearsky,firstcall,lastcall)
-
 
       use radinc_h
@@ -201,5 +200,4 @@
       giaer(:,:,:) =0.0
       radfixed=.false.
-
 
       if(firstcall) then
@@ -275,6 +273,6 @@
          endif
 
-         OLR_nu=0.
-         OSR_nu=0.
+         OLR_nu(:,:) = 0.
+         OSR_nu(:,:) = 0.
 
          if (ngridmx.eq.1) then
@@ -429,5 +427,4 @@
            reffrad,QREFvis3d,QREFir3d,                             & 
            tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths
-
 
 !-----------------------------------------------------------------------
@@ -933,3 +930,4 @@
       endif
 
+
     end subroutine callcorrk
Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 716)
@@ -8,5 +8,5 @@
      &   , iaervar,iddist,topdustref,callstats,calleofdump              &
      &   , enertest                                                     &
-     &   , callgasvis,Continuum, graybody                               &
+     &   , callgasvis,Continuum,H2Ocont_simple,graybody                 &
      &   , Nmix_co2, Nmix_h2o                                           &
      &   , dusttau                                                      &
@@ -52,5 +52,5 @@
      &   , season,diurnal,tlocked,lwrite                                &
      &   , callstats,calleofdump                                        &
-     &   , callgasvis,Continuum,graybody
+     &   , callgasvis,Continuum,H2Ocont_simple,graybody
 
       logical enertest
Index: trunk/LMDZ.GENERIC/libf/phystd/cp_neutral.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/cp_neutral.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/cp_neutral.F90	(revision 716)
@@ -8,4 +8,8 @@
   double precision T
 
+
+  ! this function has been disabled in gradients_kcm.F90 because it dont
+  ! work if you have gaseous mixtures. need to decide whether to generalise
+  ! it or simply remove entirely...
 
   ! Cp_n : cf CO2 dans abe&matsui (1988)
Index: trunk/LMDZ.GENERIC/libf/phystd/datafile_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/datafile_mod.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/datafile_mod.F90	(revision 716)
@@ -4,8 +4,10 @@
       implicit none
 
+      ! Default for Berserker @ UChicago:
+!      character(len=300) :: datadir='/home/rwordsworth/datagcm'
       ! Default for Gnome Idataplex:
-      character(len=300) :: datadir='/san/home/rdword/gcm/datagcm'
+!      character(len=300) :: datadir='/san/home/rdword/gcm/datagcm'
       ! Default for LMD machines:
-!      character(len=300) :: datadir='/u/rwlmd/datagcm'
+      character(len=300) :: datadir='/u/rwlmd/datagcm'
 
       end module datafile_mod
Index: trunk/LMDZ.GENERIC/libf/phystd/gases_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/gases_h.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/gases_h.F90	(revision 716)
@@ -5,11 +5,11 @@
 !======================================================================C
 !
-!     GASES_H    
+!     gases_h    
 !
-!     THIS A F90-ALLOCATABLE VERSION FOR gases.h -- AS 12/2011
+!     A F90-allocatable version for gases.h -- AS 12/2011
 !
 !======================================================================C
 
-      !!! THOSE ARE SET AND ALLOCATED IN su_gases.F90
+      ! Set and allocated in su_gases.F90
       integer :: ngasmx
       integer :: vgas
@@ -17,3 +17,16 @@
       real,allocatable,DIMENSION(:) :: gfrac
 
+      ! in analogy with tracer.h ...
+      integer :: igas_H2
+      integer :: igas_He
+      integer :: igas_H2O
+      integer :: igas_CO2
+      integer :: igas_CO
+      integer :: igas_N2
+      integer :: igas_O2
+      integer :: igas_SO2
+      integer :: igas_H2S
+      integer :: igas_CH4
+      integer :: igas_NH3
+
       end module gases_h
Index: trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90	(revision 716)
@@ -14,5 +14,6 @@
 
   ! internal
-  double precision cp_n,cp_v
+  !double precision cp_n,cp_v
+  double precision cp_v
 
   double precision press, rho_plus, rho_minus, dVdT, rho_c
@@ -26,6 +27,5 @@
   double precision cp_neutral
 
-  cp_n = cp_neutral(T)
-
+  !cp_n = cp_neutral(T)
 
   select case(profil_flag)
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 716)
@@ -204,4 +204,8 @@
          write(*,*) " Continuum = ",Continuum
 
+         write(*,*) "use analytic function for H2O continuum ?"
+         H2Ocont_simple=.false. ! default value
+         call getin("H2Ocont_simple",H2Ocont_simple)
+         write(*,*) " H2Ocont_simple = ",H2Ocont_simple
  
          write(*,*) "call turbulent vertical diffusion ?"
@@ -221,5 +225,5 @@
 
          write(*,*) "call CO2 condensation ?"
-         co2cond=.true. ! default value
+         co2cond=.false. ! default value
          call getin("co2cond",co2cond)
          write(*,*) " co2cond = ",co2cond
@@ -228,6 +232,6 @@
             print*,'We need a CO2 ice tracer to condense CO2'
             call abort
-         endif	 
-	 
+         endif 
+ 
          write(*,*) "CO2 supersaturation level ?"
          co2supsat=1.0 ! default value
@@ -436,5 +440,5 @@
 
          write(*,*) "Gravitationnal sedimentation ?"
-         sedimentation=.true. ! default value
+         sedimentation=.false. ! default value
          call getin("sedimentation",sedimentation)
          write(*,*) " sedimentation = ",sedimentation
Index: trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90	(revision 716)
@@ -6,9 +6,9 @@
 !     -------
 !     Calculates the H2-H2 CIA opacity, using a lookup table from
-!     Borysow (2002)
+!     HITRAN (2011)
 !
 !     Authors
 !     -------
-!     R. Wordsworth (2009)
+!     R. Wordsworth (2011)
 !     
 !==================================================================
@@ -26,6 +26,11 @@
 
       integer nS,nT
-      parameter(nS=1649)
-      parameter(nT=14)
+      parameter(nS=2428)
+      parameter(nT=10)
+
+      double precision, parameter :: losch = 2.6867774e19
+      ! Loschmit's number (molecule cm^-3 at STP) 
+      ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2
+      ! see Richard et al. 2011, JQSRT for details
 
       double precision amagat
@@ -33,7 +38,6 @@
       double precision temp_arr(nT)
       double precision abs_arr(nS,nT)
-      double precision data_tmp(nT/2+1)
-
-      integer k
+
+      integer k,iT
       logical firstcall
 
@@ -43,66 +47,56 @@
       integer strlen,ios
 
-      amagat=(273.15/temp)*(pres/101325.0)
-
-      if(firstcall)then
+      character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
+
+      character*20 bleh
+      double precision blah, Ttemp
+      integer nres
+
+      if(temp.gt.400)then
+         print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you '
+         print*,'really want to run simulations with hydrogen at T > 400 K, contact'
+         print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
+         stop
+      endif
+
+      amagat = (273.15/temp)*(pres/101325.0)
+
+      if(firstcall)then ! called by sugas_corrk only
+         print*,'----------------------------------------------------'
+         print*,'Initialising H2-H2 continuum from HITRAN database...'
 
 !     1.1 Open the ASCII files
-
-         ! cold
-         dt_file=TRIM(datadir)//'/continuum_data/H2H2_CIA_LT.dat'
+         dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
+
          open(33,file=dt_file,form='formatted',status='old',iostat=ios)
          if (ios.ne.0) then        ! file not found
            write(*,*) 'Error from interpolateH2H2'
            write(*,*) 'Data file ',trim(dt_file),' not found.'
-           write(*,*)'Check that your path to datagcm:',trim(datadir)
-           write(*,*)' is correct. You can change it in callphys.def with:'
-           write(*,*)' datadir = /absolute/path/to/datagcm'
-           write(*,*)' Also check that there is a continuum_data/H2H2_CIA_LT.dat there.'
+           write(*,*) 'Check that your path to datagcm:',trim(datadir)
+           write(*,*) 'is correct. You can change it in callphys.def with:'
+           write(*,*) 'datadir = /absolute/path/to/datagcm'
+           write(*,*) 'Also check that the continuum data continuum_data/H2-H2_norm_2011.cia is there.'
            call abort
          else
-            do k=1,nS
-               read(33,*) data_tmp
-               wn_arr(k)=data_tmp(1)
-               abs_arr(k,1:7)=data_tmp(2:8)
+
+            do iT=1,nT
+
+               read(33,fmat1) bleh,blah,blah,nres,Ttemp
+               if(nS.ne.nres)then
+                  print*,'Resolution given in file: ',trim(dt_file)
+                  print*,'is ',nres,', which does not match nS.'
+                  print*,'Please adjust nS value in interpolateH2H2.F90'
+                  stop
+               endif
+               temp_arr(iT)=Ttemp
+
+               do k=1,nS
+                  read(33,*) wn_arr(k),abs_arr(k,it)
+               end do
+
             end do
+      
          endif
          close(33)
-
-         ! hot
-         dt_file=TRIM(datadir)//'/continuum_data/H2H2_CIA_HT.dat'
-         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
-         if (ios.ne.0) then        ! file not found
-           write(*,*) 'Error from interpolateH2H2'
-           write(*,*) 'Data file ',trim(dt_file),' not found.'
-           write(*,*)'Check that your path to datagcm:',trim(datadir)
-           write(*,*)' is correct. You can change it in callphys.def with:'
-           write(*,*)' datadir = /absolute/path/to/datagcm'
-           write(*,*)' Also check that there is a continuum_data/H2H2_CIA_HT.dat there.'
-           call abort
-         else
-            do k=1,nS
-               read(34,*) data_tmp
-               wn_arr(k)=data_tmp(1)
-               ! wn_arr is identical
-               abs_arr(k,8:14)=data_tmp(2:8)
-            end do
-         endif
-         close(34)
-
-         temp_arr(1)  = 60.
-         temp_arr(2)  = 100.
-         temp_arr(3)  = 150.
-         temp_arr(4)  = 200.
-         temp_arr(5)  = 250.
-         temp_arr(6)  = 300.
-         temp_arr(7)  = 350.
-         temp_arr(8)  = 400.
-         temp_arr(9)  = 500.
-         temp_arr(10) = 600.
-         temp_arr(11) = 700.
-         temp_arr(12) = 800.
-         temp_arr(13) = 900.
-         temp_arr(14) = 1000.
-
 
          print*,'interpolateH2H2: At wavenumber ',wn,' cm^-1'
@@ -112,16 +106,36 @@
          call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
 
-         print*,'the absorption is ',abcoef,' cm^-1 amg^-2'
-
-         abcoef=abcoef*100.0*amagat**2 ! convert to m^-1
-
-         print*,'We have ',amagat,' amagats'
+         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*100.0*amagat**2 ! convert to m^-1
-         !print*,'The absorption is ',abcoef,' m^-1'
+         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
@@ -141,6 +155,6 @@
 
       integer nX,nY,i,j,a,b
-      parameter(nX=1649)
-      parameter(nY=14)
+      parameter(nX=2428)
+      parameter(nY=10)
 
       real*8 x_in,y_in,x,y,x1,x2,y1,y2
@@ -156,4 +170,5 @@
       x=x_in
       y=y_in
+
 
 !     1st check we're within the wavenumber range
@@ -180,5 +195,4 @@
          write(*,*) 'Warning from bilinearH2H2:'
          write(*,*) 'Outside continuum temperature range!'
-         write(*,*) y,y_arr(1),y_arr(nY)
          if(y.lt.y_arr(1))then
             y=y_arr(1)+0.01
@@ -187,6 +201,5 @@
             y=y_arr(nY)-0.01
          endif
-      endif
-      !else
+      else
 
 !     in the y (temperature) direction 2nd
@@ -202,5 +215,5 @@
             goto 20
          endif
-      !endif
+      endif
       
       f11=f2d_arr(a,b)
@@ -208,12 +221,7 @@
       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)
-      
+
+      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 716)
+++ trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90	(revision 716)
@@ -0,0 +1,212 @@
+     subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall)
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Calculates the H2-He CIA opacity, using a lookup table from
+!     HITRAN (2011)
+!
+!     Authors
+!     -------
+!     R. Wordsworth (2011)
+!     
+!==================================================================
+
+      use datafile_mod, only: datadir
+      implicit none
+
+      ! input
+      double precision wn                 ! wavenumber             (cm^-1)
+      double precision temp               ! temperature            (Kelvin)
+      double precision presH2             ! H2 partial pressure    (Pascals)
+      double precision presHe             ! He partial pressure    (Pascals)
+
+      ! output
+      double precision abcoef             ! absorption coefficient (m^-1)
+
+      integer nS,nT
+      parameter(nS=2428)
+      parameter(nT=10)
+
+      double precision, parameter :: losch = 2.6867774e19
+      ! Loschmit's number (molecule cm^-3 at STP) 
+      ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2
+      ! see Richard et al. 2011, JQSRT for details
+
+      double precision amagatH2
+      double precision amagatHe
+      double precision wn_arr(nS)
+      double precision temp_arr(nT)
+      double precision abs_arr(nS,nT)
+
+      integer k,iT
+      logical firstcall
+
+      save wn_arr, temp_arr, abs_arr
+
+      character*100 dt_file
+      integer strlen,ios
+
+      character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
+
+      character*20 bleh
+      double precision blah, Ttemp
+      integer nres
+
+
+      if(temp.gt.400)then
+         print*,'Your temperatures are too high for this H2-He CIA dataset.'
+         print*,'Please run mixed H2-He atmospheres below T = 400 K.'      
+         stop
+      endif
+
+      amagatH2 = (273.15/temp)*(presH2/101325.0)
+      amagatHe = (273.15/temp)*(presHe/101325.0)
+
+      if(firstcall)then ! called by sugas_corrk only
+         print*,'----------------------------------------------------'
+         print*,'Initialising H2-He continuum from HITRAN database...'
+
+!     1.1 Open the ASCII files
+         dt_file=TRIM(datadir)//'/continuum_data/H2-He_norm_2011.cia'
+
+         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+           write(*,*) 'Error from interpolateH2He'
+           write(*,*) 'Data file ',trim(dt_file),' not found.'
+           write(*,*) 'Check that your path to datagcm:',trim(datadir)
+           write(*,*) 'is correct. You can change it in callphys.def with:'
+           write(*,*) 'datadir = /absolute/path/to/datagcm'
+           write(*,*) 'Also check that the continuum data continuum_data/H2-He_norm_2011.cia is there.'
+           call abort
+         else
+
+            do iT=1,nT
+
+               read(33,fmat1) bleh,blah,blah,nres,Ttemp
+               if(nS.ne.nres)then
+                  print*,'Resolution given in file: ',trim(dt_file)
+                  print*,'is ',nres,', which does not match nS.'
+                  print*,'Please adjust nS value in interpolateH2He.F90'
+                  stop
+               endif
+               temp_arr(iT)=Ttemp
+
+               do k=1,nS
+                  read(33,*) wn_arr(k),abs_arr(k,it)
+               end do
+
+            end do
+      
+         endif
+         close(33)
+
+         print*,'interpolateH2He: At wavenumber ',wn,' cm^-1'
+         print*,'   temperature                 ',temp,' K'
+         print*,'   H2 partial pressure         ',presH2,' Pa'
+         print*,'   and He partial pressure     ',presHe,' Pa'
+
+         call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
+
+         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
+
+         ! 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/interpolateH2Ocont.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90	(revision 715)
+++ 	(revision )
@@ -1,261 +1,0 @@
-     subroutine interpolateH2Ocont(wn,temp,presS,presF,abcoef,firstcall)
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Calculates the H2O continuum opacity, using a lookup table from
-!     Clough (2005)
-!
-!     Authors
-!     -------
-!     R. Wordsworth (2011)
-!     
-!==================================================================
-
-      use datafile_mod, only: datadir
-      implicit none
-
-      ! input
-      double precision wn                 ! wavenumber             (cm^-1)
-      double precision temp               ! temperature            (Kelvin)
-      double precision presS              ! self-pressure          (Pascals)
-      double precision presF              ! foreign (air) pressure (Pascals)
-
-      ! output
-      double precision abcoef             ! absorption coefficient (m^-1)
-
-      integer nS,nT
-      parameter(nS=1001)
-      parameter(nT=11)
-
-      double precision kB
-      parameter(kB=1.3806488e-23)
-
-      double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec
-      double precision wn_arr(nS)
-      double precision temp_arr(nT)
-      double precision abs_arrS(nS,nT)
-      double precision abs_arrF(nS,nT)
-      double precision data_tmp(nT)
-
-      integer k
-      logical firstcall
-
-      save wn_arr, temp_arr, abs_arrS, abs_arrF
-
-      character*100 dt_file
-      integer strlen,ios
-
-      amagatS=(273.15/temp)*(presS/101325.0)
-      amagatF=(273.15/temp)*(presF/101325.0)
-
-      if(firstcall)then ! called by sugas_corrk only
-         print*,'----------------------------------------------------'
-         print*,'Initialising H2O continuum from MT_CKD data...'
-
-!     1.1 Open the ASCII files
-
-         ! nu array
-         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_NU.dat'
-         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
-         if (ios.ne.0) then        ! file not found
-           write(*,*) 'Error from interpolateH2O_cont'
-           write(*,*) 'Data file ',trim(dt_file),' not found.'
-           write(*,*)'Check that your path to datagcm:',trim(datadir)
-           write(*,*)' is correct. You can change it in callphys.def with:'
-           write(*,*)' datadir = /absolute/path/to/datagcm'
-           write(*,*)' Also check that there is a continuum_data/H2O_CONT_NU.dat there.'
-           call abort
-         else
-            do k=1,nS
-               read(33,*) wn_arr(k)
-            enddo
-         endif
-         close(33)
-
-         ! self broadening
-         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_SELF.dat'
-         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
-         if (ios.ne.0) then        ! file not found
-           write(*,*) 'Error from interpolateH2O_cont'
-           write(*,*) 'Data file ',trim(dt_file),' not found.'
-           write(*,*)'Check that your path to datagcm:',trim(datadir)
-           write(*,*)' is correct. You can change it in callphys.def with:'
-           write(*,*)' datadir = /absolute/path/to/datagcm'
-           write(*,*)' Also check that there is a continuum_data/H2O_CONT_SELF.dat there.'
-           call abort
-         else
-            do k=1,nS
-               read(34,*) data_tmp
-               abs_arrS(k,1:nT)=data_tmp(1:nT)
-            end do
-         endif
-         close(34)
-
-         ! foreign (N2+O2+Ar) broadening
-         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_FOREIGN.dat'
-         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
-         if (ios.ne.0) then        ! file not found
-           write(*,*) 'Error from interpolateH2O_cont'
-           write(*,*) 'Data file ',trim(dt_file),' not found.'
-           write(*,*)'Check that your path to datagcm:',trim(datadir)
-           write(*,*)' is correct. You can change it in callphys.def with:'
-           write(*,*)' datadir = /absolute/path/to/datagcm'
-           write(*,*)' Also check that there is a continuum_data/H2O_CONT_FOREIGN.dat there.'
-           call abort
-         else
-            do k=1,nS
-               read(35,*) data_tmp
-               abs_arrF(k,1:nT)=data_tmp(1:nT)
-            end do
-         endif
-         close(35)
-
-         temp_arr(1)  = 200.
-         temp_arr(2)  = 250.
-         temp_arr(3)  = 300.
-         temp_arr(4)  = 350.
-         temp_arr(5)  = 400.
-         temp_arr(6)  = 450.
-         temp_arr(7)  = 500.
-         temp_arr(8)  = 550.
-         temp_arr(9)  = 600.
-         temp_arr(10) = 650.
-         temp_arr(11) = 700.
-
-         print*,'interpolateH2Ocont: At wavenumber ',wn,' cm^-1'
-         print*,'   temperature ',temp,' K'
-         print*,'   H2O pressure ',presS,' Pa'
-         print*,'   air pressure ',presF,' Pa'
-
-         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
-         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'
-
-         print*,'We have ',amagatS,' amagats of H2O vapour'
-         print*,'and ',amagatF,' amagats of air'
-
-         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
-
-         Nmolec = (presS+presF)/(kB*temp)
-
-         print*,'Total number of molecules per m^3 is',Nmolec
-
-         abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
-         abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account
-
-         print*,'So the total absorption is ',abcoef,' m^-1'
-         print*,'And optical depth / km : ',1000.0*abcoef
-
-      else
-
-         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
-         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
-
-         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
-
-         !abcoef = (presS/100000.0)*3.0e-24
-         !print*,'Matsui TEST'
-
-         Nmolec = (presS+presF)/(kB*temp)  ! sure this is correct??
-         abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
-         abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account
-
-
-         ! 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
-
-
-
-
-!-------------------------------------------------------------------------
-      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/interpolateH2Ocont_CKD.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_CKD.F90	(revision 716)
+++ trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_CKD.F90	(revision 716)
@@ -0,0 +1,259 @@
+     subroutine interpolateH2Ocont_CKD(wn,temp,presS,presF,abcoef,firstcall)
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Calculates the H2O continuum opacity, using a lookup table from
+!     Clough (2005).
+!
+!     Authors
+!     -------
+!     R. Wordsworth (2011)
+!     
+!==================================================================
+
+      use datafile_mod, only: datadir
+      implicit none
+
+      ! input
+      double precision wn                 ! wavenumber             (cm^-1)
+      double precision temp               ! temperature            (Kelvin)
+      double precision presS              ! self-pressure          (Pascals)
+      double precision presF              ! foreign (air) pressure (Pascals)
+
+      ! output
+      double precision abcoef             ! absorption coefficient (m^-1)
+
+      integer nS,nT
+      parameter(nS=1001)
+      parameter(nT=11)
+
+      double precision kB
+      parameter(kB=1.3806488e-23)
+
+      double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec
+      double precision wn_arr(nS)
+      double precision temp_arr(nT)
+      double precision abs_arrS(nS,nT)
+      double precision abs_arrF(nS,nT)
+      double precision data_tmp(nT)
+
+      integer k
+      logical firstcall
+
+      save wn_arr, temp_arr, abs_arrS, abs_arrF
+
+      character*100 dt_file
+      integer strlen,ios
+
+      amagatS=(273.15/temp)*(presS/101325.0)
+      amagatF=(273.15/temp)*(presF/101325.0)
+
+      if(firstcall)then ! called by sugas_corrk only
+         print*,'----------------------------------------------------'
+         print*,'Initialising H2O continuum from MT_CKD data...'
+
+!     1.1 Open the ASCII files
+
+         ! nu array
+         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_NU.dat'
+         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+           write(*,*) 'Error from interpolateH2O_cont'
+           write(*,*) 'Data file ',trim(dt_file),' not found.'
+           write(*,*)'Check that your path to datagcm:',trim(datadir)
+           write(*,*)' is correct. You can change it in callphys.def with:'
+           write(*,*)' datadir = /absolute/path/to/datagcm'
+           write(*,*)' Also check that there is a continuum_data/H2O_CONT_NU.dat there.'
+           call abort
+         else
+            do k=1,nS
+               read(33,*) wn_arr(k)
+            enddo
+         endif
+         close(33)
+
+         ! self broadening
+         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_SELF.dat'
+         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+           write(*,*) 'Error from interpolateH2O_cont'
+           write(*,*) 'Data file ',trim(dt_file),' not found.'
+           write(*,*)'Check that your path to datagcm:',trim(datadir)
+           write(*,*)' is correct. You can change it in callphys.def with:'
+           write(*,*)' datadir = /absolute/path/to/datagcm'
+           write(*,*)' Also check that there is a continuum_data/H2O_CONT_SELF.dat there.'
+           call abort
+         else
+            do k=1,nS
+               read(34,*) data_tmp
+               abs_arrS(k,1:nT)=data_tmp(1:nT)
+            end do
+         endif
+         close(34)
+
+         ! foreign (N2+O2+Ar) broadening
+         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_FOREIGN.dat'
+         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
+         if (ios.ne.0) then        ! file not found
+           write(*,*) 'Error from interpolateH2O_cont'
+           write(*,*) 'Data file ',trim(dt_file),' not found.'
+           write(*,*)'Check that your path to datagcm:',trim(datadir)
+           write(*,*)' is correct. You can change it in callphys.def with:'
+           write(*,*)' datadir = /absolute/path/to/datagcm'
+           write(*,*)' Also check that there is a continuum_data/H2O_CONT_FOREIGN.dat there.'
+           call abort
+         else
+            do k=1,nS
+               read(35,*) data_tmp
+               abs_arrF(k,1:nT)=data_tmp(1:nT)
+            end do
+         endif
+         close(35)
+
+         temp_arr(1)  = 200.
+         temp_arr(2)  = 250.
+         temp_arr(3)  = 300.
+         temp_arr(4)  = 350.
+         temp_arr(5)  = 400.
+         temp_arr(6)  = 450.
+         temp_arr(7)  = 500.
+         temp_arr(8)  = 550.
+         temp_arr(9)  = 600.
+         temp_arr(10) = 650.
+         temp_arr(11) = 700.
+
+         print*,'interpolateH2Ocont: At wavenumber ',wn,' cm^-1'
+         print*,'   temperature ',temp,' K'
+         print*,'   H2O pressure ',presS,' Pa'
+         print*,'   air pressure ',presF,' Pa'
+
+         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
+         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'
+
+         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
+                                                    ! 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
+
+         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
+
+      else
+
+         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
+         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
+
+         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/interpolateH2Ocont_PPC.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_PPC.F90	(revision 716)
+++ trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_PPC.F90	(revision 716)
@@ -0,0 +1,84 @@
+     subroutine interpolateH2Ocont_PPC(wn,temp,presS,presF,abcoef,firstcall)
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Calculates the H2O continuum opacity, using the formulae
+!     provided in Pierrehumbert, PPC (2010). As this is based on
+!     the CKD continuum, it provides a useful check for the 
+!     implementation of the more general interpolateH2Ocont_CKD.F90.
+!
+!     Authors
+!     -------
+!     R. Wordsworth (2012)
+!     
+!==================================================================
+
+      use watercommon_h, only: mH2O
+      use datafile_mod, only: datadir
+      implicit none
+
+      ! input
+      double precision wn                 ! wavenumber             (cm^-1)
+      double precision temp               ! temperature            (Kelvin)
+      double precision presS              ! self-pressure          (Pascals)
+      double precision presF              ! foreign (air) pressure (Pascals)
+
+      ! parameters
+      double precision, parameter :: T0 = 296.0
+      double precision, parameter :: p0 = 1.D+4
+
+      ! variables
+      double precision rho_w, x
+
+      ! output
+      double precision abcoef             ! absorption coefficient (m^-1)
+
+      logical firstcall
+
+      x = wn - 2500.
+
+      if(firstcall)then ! called by sugas_corrk only
+         print*,'----------------------------------------------------'
+         print*,'Testing H2O continuum...'
+
+         print*,'interpolateH2Ocont: At wavenumber ',wn,' cm^-1'
+         print*,'   temperature ',temp,' K'
+         print*,'   H2O pressure ',presS,' Pa'
+
+         rho_w = presS/((8.31446/(mH2O/1000.))*temp)
+
+         if(wn.gt.500 .and. wn.lt.1400)then
+            abcoef = exp(12.167 - 0.050898*wn + 8.3207e-5*wn**2 - 7.0748e-8*wn**3 + 2.3261e-11*wn**4)*(T0/temp)**4.25*(presS/p0)
+         elseif(wn.gt.2100 .and. wn.lt.3000)then
+            abcoef = exp(-6.0055 - 0.0021363*x + 6.4723e-7*x**2 - 1.493e-8*x**3 + 2.5621e-11*x**4 + 7.328e-14*x**5)*(T0/temp)**4.25*(presS/p0)
+         else
+            abcoef = 0.0
+         endif
+         abcoef = abcoef*rho_w
+
+         print*,'The self absorption is ',abcoef,' m^-1'
+         print*,'And optical depth / km : ',1000.0*abcoef
+
+      else
+
+         rho_w = presS/((8.31446/(mH2O/1000.))*temp)
+
+         if(wn.gt.500 .and. wn.lt.1400)then
+            abcoef = exp(12.167 - 0.050898*wn + 8.3207e-5*wn**2 - 7.0748e-8*wn**3 + 2.3261e-11*wn**4)*(T0/temp)**4.25*(presS/p0)
+         elseif(wn.gt.2100 .and. wn.lt.3000)then
+            abcoef = exp(-6.0055 - 0.0021363*x + 6.4723e-7*x**2 - 1.493e-8*x**3 + 2.5621e-11*x**4 + 7.328e-14*x**5)*(T0/temp)**4.25*(presS/p0)
+         else
+            abcoef = 0.0
+         endif
+         abcoef = abcoef*rho_w
+
+         ! 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_PPC
+
Index: trunk/LMDZ.GENERIC/libf/phystd/interpolateN2H2.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/interpolateN2H2.F90	(revision 716)
+++ trunk/LMDZ.GENERIC/libf/phystd/interpolateN2H2.F90	(revision 716)
@@ -0,0 +1,209 @@
+subroutine interpolateN2H2(wn,temp,presN2,presH2,abcoef,firstcall)
+
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Calculates the N2-H2 CIA opacity, using a lookup table from
+  !     HITRAN (2011)
+  !
+  !     Authors
+  !     -------
+  !     R. Wordsworth (2011)
+  !     
+  !==================================================================
+
+  use datafile_mod, only: datadir
+  implicit none
+
+  ! input
+  double precision wn                 ! wavenumber             (cm^-1)
+  double precision temp               ! temperature            (Kelvin)
+  double precision presN2             ! N2 partial pressure    (Pascals)
+  double precision presH2             ! H2 partial pressure    (Pascals)
+
+  ! output
+  double precision abcoef             ! absorption coefficient (m^-1)
+
+  integer nS,nT
+  parameter(nS=1914)
+  parameter(nT=10)
+
+  double precision, parameter :: losch = 2.6867774e19
+  ! Loschmit's number (molecule cm^-3 at STP) 
+  ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2
+  ! see Richard et al. 2011, JQSRT for details
+
+  double precision amagatN2
+  double precision amagatH2
+  double precision wn_arr(nS)
+  double precision temp_arr(nT)
+  double precision abs_arr(nS,nT)
+
+  integer k,iT
+  logical firstcall
+
+  save wn_arr, temp_arr, abs_arr
+
+  character*100 dt_file
+  integer strlen,ios
+
+  character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
+
+  character*20 bleh
+  double precision blah, Ttemp
+  integer nres
+
+  if(temp.gt.400)then
+     print*,'Your temperatures are too high for this N2-H2 CIA dataset.'
+     print*,'Please run mixed N2-H2 atmospheres below T = 400 K.'      
+     stop
+  endif
+
+  amagatN2 = (273.15/temp)*(presN2/101325.0)
+  amagatH2 = (273.15/temp)*(presH2/101325.0)
+
+  if(firstcall)then ! called by sugas_corrk only
+     print*,'----------------------------------------------------'
+     print*,'Initialising N2-H2 continuum from HITRAN database...'
+
+     !     1.1 Open the ASCII files
+     dt_file=TRIM(datadir)//'/continuum_data/N2-H2_2011.cia'
+
+     open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+     if (ios.ne.0) then        ! file not found
+        write(*,*) 'Error from interpolateN2H2'
+        write(*,*) 'Data file ',trim(dt_file),' not found.'
+        write(*,*) 'Check that your path to datagcm:',trim(datadir)
+        write(*,*) 'is correct. You can change it in callphys.def with:'
+        write(*,*) 'datadir = /absolute/path/to/datagcm'
+        write(*,*) 'Also check that the continuum data continuum_data/N2-H2_2011.cia is there.'
+        call abort
+     else
+
+        do iT=1,nT
+
+           read(33,fmat1) bleh,blah,blah,nres,Ttemp
+           if(nS.ne.nres)then
+              print*,'Resolution given in file: ',trim(dt_file)
+              print*,'is ',nres,', which does not match nS.'
+              print*,'Please adjust nS value in interpolateN2H2.F90'
+              stop
+           endif
+           temp_arr(iT)=Ttemp
+
+           do k=1,nS
+              read(33,*) wn_arr(k),abs_arr(k,it)
+           end do
+
+        end do
+
+     endif
+     close(33)
+
+     print*,'interpolateN2H2: At wavenumber ',wn,' cm^-1'
+     print*,'   temperature                 ',temp,' K'
+     print*,'   N2 partial pressure         ',presN2,' Pa'
+     print*,'   and H2 partial pressure     ',presH2,' Pa'
+
+     call bilinearN2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
+
+     print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
+     print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
+
+     abcoef=abcoef*losch**2*100.0*amagatN2*amagatH2 ! convert to m^-1
+
+     print*,'We have ',amagatN2,' amagats of N2'
+     print*,'and     ',amagatH2,' amagats of H2'
+     print*,'So the absorption is ',abcoef,' m^-1'
+
+  else
+
+     call bilinearN2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
+     abcoef=abcoef*losch**2*100.0*amagatN2*amagatH2 ! 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 interpolateN2H2
+
+
+!-------------------------------------------------------------------------
+subroutine bilinearN2H2(x_arr,y_arr,f2d_arr,x_in,y_in,f)
+
+  implicit none
+
+  integer nX,nY,i,j,a,b
+  parameter(nX=1914)
+  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 bilinearN2H2:'
+     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 bilinearN2H2
Index: trunk/LMDZ.GENERIC/libf/phystd/interpolateN2N2.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/interpolateN2N2.F90	(revision 716)
+++ trunk/LMDZ.GENERIC/libf/phystd/interpolateN2N2.F90	(revision 716)
@@ -0,0 +1,205 @@
+subroutine interpolateN2N2(wn,temp,pres,abcoef,firstcall)
+
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Calculates the N2-N2 CIA opacity, using a lookup table from
+  !     HITRAN (2011)
+  !
+  !     Authors
+  !     -------
+  !     R. Wordsworth (2011)
+  !     
+  !==================================================================
+
+  use datafile_mod, only: datadir
+  implicit none
+
+  ! input
+  double precision wn                 ! wavenumber             (cm^-1)
+  double precision temp               ! temperature            (Kelvin)
+  double precision pres               ! pressure               (Pascals)
+
+  ! output
+  double precision abcoef             ! absorption coefficient (m^-1)
+
+  integer nS,nT
+  parameter(nS=582)
+  parameter(nT=10)
+
+  double precision, parameter :: losch = 2.6867774e19
+  ! Loschmit's number (molecule cm^-3 at STP) 
+  ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2
+  ! see Richard et al. 2011, JQSRT for details
+
+  double precision amagat
+  double precision wn_arr(nS)
+  double precision temp_arr(nT)
+  double precision abs_arr(nS,nT)
+
+  integer k,iT
+  logical firstcall
+
+  save wn_arr, temp_arr, abs_arr
+
+  character*100 dt_file
+  integer strlen,ios
+
+  character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
+
+  character*20 bleh
+  double precision blah, Ttemp
+  integer nres
+
+
+  if(temp.gt.400)then
+     print*,'Your temperatures are too high for this N2-N2 CIA dataset.'
+     print*,'Currently, HITRAN provides data for this pair in the range 40-400 K.'
+     stop
+  endif
+
+  amagat = (273.15/temp)*(pres/101325.0)
+
+  if(firstcall)then ! called by sugas_corrk only
+     print*,'----------------------------------------------------'
+     print*,'Initialising N2-N2 continuum from HITRAN database...'
+
+     !     1.1 Open the ASCII files
+     dt_file=TRIM(datadir)//'/continuum_data/N2-N2_2011.cia'
+
+     open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+     if (ios.ne.0) then        ! file not found
+        write(*,*) 'Error from interpolateN2N2'
+        write(*,*) 'Data file ',trim(dt_file),' not found.'
+        write(*,*) 'Check that your path to datagcm:',trim(datadir)
+        write(*,*) 'is correct. You can change it in callphys.def with:'
+        write(*,*) 'datadir = /absolute/path/to/datagcm'
+        write(*,*) 'Also check that the continuum data continuum_data/N2-N2_norm_2011.cia is there.'
+        call abort
+     else
+
+        do iT=1,nT
+
+           read(33,fmat1) bleh,blah,blah,nres,Ttemp
+           if(nS.ne.nres)then
+              print*,'Resolution given in file: ',trim(dt_file)
+              print*,'is ',nres,', which does not match nS.'
+              print*,'Please adjust nS value in interpolateN2N2.F90'
+              stop
+           endif
+           temp_arr(iT)=Ttemp
+
+           do k=1,nS
+              read(33,*) wn_arr(k),abs_arr(k,it)
+           end do
+
+        end do
+
+     endif
+     close(33)
+
+     print*,'interpolateN2N2: At wavenumber ',wn,' cm^-1'
+     print*,'   temperature ',temp,' K'
+     print*,'   pressure ',pres,' Pa'
+
+     call bilinearN2N2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
+
+     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 N2'
+     print*,'So the absorption is ',abcoef,' m^-1'
+
+  else
+
+     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/kcm1d.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 716)
@@ -8,21 +8,21 @@
   implicit none
 
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Run the universal model radiative transfer once in a 1D column.
-!     Useful for climate evolution studies etc.
-!     
-!     It can be compiled with a command like (e.g. for 25 layers):
-!                                  "makegcm -p std -d 25 kcm1d"
-!     It requires the files "callphys.def", "gases.def"
-!     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
-!
-!     Authors
-!     -------
-!     R. Wordsworth
-!
-!==================================================================
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Run the universal model radiative transfer once in a 1D column.
+  !     Useful for climate evolution studies etc.
+  !     
+  !     It can be compiled with a command like (e.g. for 25 layers):
+  !                                  "makegcm -p std -d 25 kcm1d"
+  !     It requires the files "callphys.def", "gases.def"
+  !     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
+  !
+  !     Authors
+  !     -------
+  !     R. Wordsworth
+  !
+  !==================================================================
 
 #include "dimensions.h"
@@ -35,7 +35,7 @@
 #include "advtrac.h"
 
-! --------------------------------------------------------------
-!  Declarations
-! --------------------------------------------------------------
+  ! --------------------------------------------------------------
+  !  Declarations
+  ! --------------------------------------------------------------
 
   integer nlayer,nlevel,nq
@@ -61,6 +61,6 @@
   real fluxabs_sw    ! SW flux absorbed by planet (W/m2)
   real fluxtop_dn    ! incident top of atmosphere SW flux (W/m2)
-      
-! not used
+
+  ! not used
   real reffrad(nlayermx,naerkind)
   real nueffrad(nlayermx,naerkind)
@@ -71,14 +71,14 @@
   real dTstrat
   real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg)
-  real OLR_nu(L_NSPECTI)
-  real GSR_nu(L_NSPECTV)
+  real OLR_nu(ngridmx,L_NSPECTI)
+  real OSR_nu(ngridmx,L_NSPECTV)
   real Eatmtot
 
   integer ierr
-  logical firstcall,lastcall
-
-! --------------
-! Initialisation
-! --------------
+  logical firstcall,lastcall,global1d
+
+  ! --------------
+  ! Initialisation
+  ! --------------
 
   pi=2.E+0*asin(1.E+0)
@@ -89,13 +89,18 @@
   totcloudfrac  = 0.0
 
-!  Load parameters from "run.def" 
-! -------------------------------
-
-! check if 'run1d.def' file is around (otherwise reading parameters
-! from callphys.def via getin() routine won't work.)
-  open(90,file='run.def',status='old',form='formatted',&
+
+  nlayer=nlayermx
+  nlevel=nlayer+1
+
+
+  !  Load parameters from "run.def" 
+  ! -------------------------------
+
+  ! check if 'kcm1d.def' file is around (otherwise reading parameters
+  ! from callphys.def via getin() routine won't work.)
+  open(90,file='kcm1d.def',status='old',form='formatted',&
        iostat=ierr)
   if (ierr.ne.0) then
-     write(*,*) 'Cannot find required file "run.def"'
+     write(*,*) 'Cannot find required file "kcm1d.def"'
      write(*,*) '  (which should contain some input parameters'
      write(*,*) '   along with the following line:'
@@ -108,6 +113,28 @@
   endif
 
-  nlayer=nlayermx
-  nlevel=nlayer+1
+! now, run.def is needed anyway. so we create a dummy temporary one
+! for ioipsl to work. if a run.def is already here, stop the
+! program and ask the user to do a bit of cleaning
+  open(90,file='run.def',status='old',form='formatted',& 
+       iostat=ierr)
+  if (ierr.eq.0) then
+     close(90)
+     write(*,*) 'There is already a run.def file.'
+     write(*,*) 'This is not compatible with 1D runs.'
+     write(*,*) 'Please remove the file and restart the run.'
+     write(*,*) 'Runtime parameters are supposed to be in kcm1d.def'
+     stop
+  else
+     call system('touch run.def')
+     call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
+     call system("echo 'INCLUDEDEF=kcm1d.def' >> run.def")
+  endif
+
+  global1d = .false. ! default value
+  call getin("global1d",global1d)
+  if(.not.global1d)then
+     print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
+     stop
+  end if
 
   psurf_n=100000. ! default value for psurf
@@ -115,5 +142,9 @@
   call getin("psurf",psurf_n)
   write(*,*) " psurf = ",psurf_n
-      
+
+! OK. now that run.def has been read once -- any variable is in memory.
+! so we can dump the dummy run.def
+  call system("rm -rf run.def")
+
   tsurf=300.0 ! default value for tsurf
   print*,'Surface temperature (K)?'
@@ -125,5 +156,5 @@
   call getin("g",g)
   write(*,*) " g = ",g
-  
+
   periastr = 1.0
   apoastr  = 1.0
@@ -136,5 +167,5 @@
   call getin("apoastr",apoastr)
   write(*,*) "apoastron = ",apoastr
- 
+
   albedo=0.2 ! default value for albedo
   print*,'Albedo of bare ground?'
@@ -155,12 +186,19 @@
   cpp=0.
 
+  check_cpp_match = .false.
+  call getin("check_cpp_match",check_cpp_match)
+  if (check_cpp_match) then
+     print*,"In 1D modeling, check_cpp_match is supposed to be F"
+     print*,"Please correct callphys.def"
+     stop
+  endif
+
   call su_gases
   call calc_cpp_mugaz
 
-
   call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp)
 
-! Tracer initialisation
-! ---------------------
+  ! Tracer initialisation
+  ! ---------------------
   if (tracer) then
      ! load tracer names from file 'traceur.def'
@@ -183,5 +221,5 @@
            endif
         endif
-        
+
         do iq=1,nq
            ! minimal version, just read in the tracer names, 1 per line
@@ -193,7 +231,6 @@
         enddo !of do iq=1,nq
      endif
-  !endif
-  
-  call initracer()
+
+     call initracer()
 
   endif
@@ -205,5 +242,5 @@
      enddo
   enddo
-  
+
   do iq=1,nqmx
      qsurf(iq) = 0.
@@ -214,26 +251,26 @@
 
   iter    = 1
-  Tstrat  = 60.0
+  Tstrat  = 150.0
   dTstrat = 1000.0
 
-! ---------
-! Run model
-! ---------
-  do
+  ! ---------
+  ! Run model
+  ! ---------
+  !do
      psurf = psurf_n
 
-!    Create vertical profiles
+     !    Create vertical profiles
      call kcmprof_fn(psurf,qsurf(1),tsurf,    &
           tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
 
-!    Run radiative transfer
+     !    Run radiative transfer
      call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
           albedo,emis,mu0,plev,play,temp,                    &
-          tsurf,fract,dist_star,aerosol,cpp,muvar,         &
+          tsurf,fract,dist_star,aerosol,muvar,         &
           dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-          fluxabs_sw,fluxtop_dn,reffrad,tau_col,  &
+          fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,reffrad,tau_col,  &
           cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
- 
-!    Iterate stratospheric temperature
+
+     !    Iterate stratospheric temperature
      print*,'Tstrat = ',Tstrat
      dTstrat = Tstrat
@@ -244,26 +281,38 @@
      dTstrat = dTstrat-Tstrat
 
-     if(abs(dTstrat).lt.1.0)then
-        print*,'dTstrat = ',dTstrat
-        exit
-     endif
-
-     iter=iter+1
-     if(iter.eq.100)then
-        print*,'Stratosphere failed to converge after'
-        print*,'100 iteration, aborting run.'
-        call abort
-     endif
-
-  end do
-
-! Calculate total atmospheric energy
+     !if(abs(dTstrat).lt.1.0)then
+     !   print*,'dTstrat = ',dTstrat
+     !   exit
+     !endif
+
+     !iter=iter+1
+     !if(iter.eq.100)then
+     !   print*,'Stratosphere failed to converge after'
+     !   print*,'100 iteration, aborting run.'
+     !   call abort
+     !endif
+
+  !end do
+
+  ! Run radiative transfer one last time to get OLR,OSR
+  firstcall=.false.
+  lastcall=.true.
+  call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
+       albedo,emis,mu0,plev,play,temp,                    &
+       tsurf,fract,dist_star,aerosol,muvar,         &
+       dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
+       fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,             &
+       reffrad,tau_col,  &
+       cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
+
+
+  ! Calculate total atmospheric energy
   Eatmtot=0.0
-!  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
-!     q(:,1),muvar,Eatmtot)
-
-! ------------------------
-! Save data to ascii files
-! ------------------------
+  !  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
+  !     q(:,1),muvar,Eatmtot)
+
+  ! ------------------------
+  ! Save data to ascii files
+  ! ------------------------
 
   print*,'Saving profiles...'
@@ -293,4 +342,6 @@
   close(119)
 
+  print*, tsurf,psurf,fluxtop_dn,fluxabs_sw,fluxtop_lw 
+
   print*,'Saving scalars...'
   open(116,file='surf_vals.out')
@@ -305,11 +356,11 @@
   open(117,file='OLRnu.out')
   do iw=1,L_NSPECTI
-     write(117,*) OLR_nu(iw)
+     write(117,*) OLR_nu(1,iw)
   enddo
   close(117)
-  
-  open(127,file='GSRnu.out')
+
+  open(127,file='OSRnu.out')
   do iw=1,L_NSPECTV
-     write(127,*) GSR_nu(iw)
+     write(127,*) OSR_nu(1,iw)
   enddo
   close(127)  
Index: trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 716)
@@ -40,7 +40,5 @@
   double precision Dz, Dp
   double precision Ptop, dlogp, Psat_max
-  parameter (Ptop=1.0)           ! Pressure at TOA [Pa]
-  !parameter (Psat_max=100000.0) ! Maximum vapour pressure [Pa]
-  parameter (Psat_max=0.0) ! set to zero for dry calculations
+  parameter (Ptop=1.0)                             ! Pressure at TOA [Pa]
 
   double precision T(1:nlay)                       ! temperature [K]
@@ -64,6 +62,7 @@
   double precision psat_v      ! local Psat_H2O value
   double precision Tcrit       ! Critical temperature [K]
-!  double precision Tsat        ! local Tsat_H2O value
   double precision rho_vTEMP,rho_nTEMP
+
+  double precision TCO2cond ! for CO2 condensation quasi-hack
 
   ! variables necessary for steam.f90
@@ -74,5 +73,8 @@
 
   logical verbose
-  parameter(verbose=.false.)
+  parameter(verbose=.true.)
+
+  logical add_Pvar_to_total
+  parameter(add_Pvar_to_total=.true.)
 
   ! initialise flags
@@ -82,24 +84,26 @@
   ! assign input variables
   m_n     = dble(mugaz/1000.)
+  cp_n    = cpp
   ! modify/generalise later??
 
-  if(ngasmx.gt.2)then
-     print*,'Only two species possible at present'
-     stop
-  elseif(ngasmx.eq.1)then
+  Psat_max = 1000000.0 ! maximum vapour pressure [Pa]
+                       ! set huge until further notice 
+
+  if(vgas.lt.1)then
      if(psat_max.gt.0.0)then
         print*,'Must have Psat_max=0 if no variable species'
-        stop
+        psat_max=0.0
+        !stop
      endif
      print*, 'Assuming pure atmosphere'
      m_v   = 1.0
      tcrit = 1000.0
-  elseif(gnom(ngasmx).eq.'H2O')then
+  elseif(gnom(vgas).eq.'H2O')then
      m_v   = dble(mH2O/1000.)
      tcrit = 6.47d2
-  elseif(gnom(ngasmx).eq.'NH3')then
+  elseif(gnom(vgas).eq.'NH3')then
      m_v   = 17.031/1000.
      tcrit = 4.06d2
-  elseif(gnom(ngasmx).eq.'CH4')then
+  elseif(gnom(vgas).eq.'CH4')then
      m_v   = 16.04/1000.
      tcrit = 1.91d2
@@ -111,24 +115,14 @@
 
   rmn     = rc/m_n
-  
-  Psurf_n = dble(psurf_rcm)
   Ttop    = dble(Tstra_rcm)
   Tsurf   = dble(Tsurf_rcm)
 
-
-  ! assume vapour saturation at surface (for now)
-!  if (tsurf > tcrit) then
-!     print*,'Above critical temperature for ',gnom(2),&
-!          ', at surface, assuming 10 bar pressure instead'
-!     Psurf_v = 10*1d5
-!     profil_flag(1) = 0
-!i!  endif
-
-
-  psat_v=psat_max
-  if(gnom(ngasmx).eq.'H2O')then
-     call Psat_H2O(tsurf,psat_v)
-  elseif(gnom(ngasmx).eq.'NH3')then
-     call Psat_NH3(tsurf,psat_v)
+  psat_v  = psat_max
+  if(vgas.gt.0)then
+     if(gnom(vgas).eq.'H2O')then
+        call Psat_H2O(tsurf,psat_v)
+     elseif(gnom(vgas).eq.'NH3')then
+        call Psat_NH3(tsurf,psat_v)
+     endif
   endif
 
@@ -142,5 +136,16 @@
   endif
 
-
+  if(add_Pvar_to_total)then
+    Psurf_n = dble(psurf_rcm)
+    psurf_rcm = real(Psurf_n+Psurf_v)
+  else
+    Psurf_n = dble(psurf_rcm) - Psurf_v
+  endif
+
+  ! include relative humidity option
+  !if(satval.lt.1.0)then
+  !   Psurf_v = Psurf_v*satval
+  !   profil_flag(1) = 0     
+  !endif
 
   if(verbose)then
@@ -155,8 +160,5 @@
   endif
 
-
-
   ! define fine pressure grid
-  psurf_rcm = real(Psurf_n+Psurf_v)
   dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayermx
 
@@ -207,7 +209,4 @@
   do ilay=2,nlay-1
 
-
-
-
      ! calculate altitude levels (for diagnostic only)
      Dz         = -Dp/(  g*(rho_n(ilay) + rho_v(ilay))  )
@@ -228,8 +227,11 @@
      ! test for moist adiabat at next level
      psat_v=psat_max
-     if(gnom(ngasmx).eq.'H2O')then
+
+     if(vgas.gt.0)then
+     if(gnom(vgas).eq.'H2O')then
         call Psat_H2O(T(ilay+1),psat_v)
-     elseif(gnom(ngasmx).eq.'NH3')then
+     elseif(gnom(vgas).eq.'NH3')then
         call Psat_NH3(T(ilay+1),psat_v)
+     endif
      endif
 
@@ -253,8 +255,11 @@
        
         psat_v=psat_max 
-        if(gnom(ngasmx).eq.'H2O')then
+
+        if(vgas.gt.0)then
+        if(gnom(vgas).eq.'H2O')then
            call Psat_H2O(T(ilay+1),psat_v)
-        elseif(gnom(ngasmx).eq.'NH3')then
+        elseif(gnom(vgas).eq.'NH3')then
            call Psat_NH3(T(ilay+1),psat_v)
+        endif
         endif
 
@@ -281,17 +286,4 @@
      end select
 
-
-
-!         print*,'profil_flag=',profil_flag(ilay)
-!         print*,'T=',T(ilay)
-!         print*,'a=',rho_v(ilay)/rho_n(ilay)
-!      if(profil_flag(ilay).eq.1)then
-!           stop
-!     endif
-
-
-
-
-
   enddo
 
@@ -378,4 +370,22 @@
   enddo
 
+!    CO2 condensation 'haircut' of temperature profile if necessary
+  if(co2cond)then
+     print*,'CO2 condensation haircut - assumes CO2-dominated atmosphere!'
+     do ilay=2,nlayermx
+        if(P_rcm(ilay).lt.518000.)then
+           TCO2cond = (-3167.8)/(log(.01*P_rcm(ilay))-23.23) ! Fanale's formula
+        else
+           TCO2cond = 684.2-92.3*log(P_rcm(ilay))+4.32*log(P_rcm(ilay))**2 
+           ! liquid-vapour transition (based on CRC handbook 2003 data)
+        endif
+
+        print*,'p=',P_rcm(ilay),', T=',T_rcm(ilay),' Tcond=',TCO2cond
+        if(T_rcm(ilay).lt.TCO2cond)then
+           T_rcm(ilay)=TCO2cond
+        endif
+     enddo
+  endif
+
   return
 end subroutine kcmprof_fn
Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 716)
@@ -1,31 +1,30 @@
-      subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
-           QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
-           TMID,PMID,TAUGSURF,QVAR,MUVAR)
-
-      use radinc_h
-      use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi,scalep
-      use gases_h
-      implicit none
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Calculates longwave optical constants at each level. For each
-!     layer and spectral interval in the IR it calculates WBAR, DTAU
-!     and COSBAR. For each level it calculates TAU.
-!     
-!     TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
-!     at the *bottom* of layer L), LW is the spectral wavelength interval.
-!     
-!     TLEV(L) - Temperature at the layer boundary (i.e., level)
-!     PLEV(L) - Pressure at the layer boundary (i.e., level)
-!
-!     Authors
-!     -------
-!     Adapted from the NASA Ames code by R. Wordsworth (2009)
-!     
-!==================================================================
-
+subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
+     QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
+     TMID,PMID,TAUGSURF,QVAR,MUVAR)
+
+  use radinc_h
+  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep
+  use gases_h
+  implicit none
+
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Calculates longwave optical constants at each level. For each
+  !     layer and spectral interval in the IR it calculates WBAR, DTAU
+  !     and COSBAR. For each level it calculates TAU.
+  !     
+  !     TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
+  !     at the *bottom* of layer L), LW is the spectral wavelength interval.
+  !     
+  !     TLEV(L) - Temperature at the layer boundary (i.e., level)
+  !     PLEV(L) - Pressure at the layer boundary (i.e., level)
+  !
+  !     Authors
+  !     -------
+  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
+  !     
+  !==================================================================
 
 #include "comcstfi.h"
@@ -33,297 +32,351 @@
 
 
-      real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
-      real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
-      real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
-      real*8 PLEV(L_LEVELS)
-      real*8 TLEV(L_LEVELS)
-      real*8 TMID(L_LEVELS), PMID(L_LEVELS)
-      real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-
-!     For aerosols
-      real*8  QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
-      real*8  QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
-      real*8  GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
-      real*8  TAUAERO(L_LEVELS+1,NAERKIND)
-      real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
-      real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
-
-      integer L, NW, NG, K, LK, IAER
-      integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
-      real*8  ANS, TAUGAS
-      real*8  DPR(L_LEVELS), U(L_LEVELS)
-      real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
-
-      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
-      real*8 DCONT
-      double precision wn_cont, p_cont, p_air, T_cont, dtemp
-
-!     Variable species mixing ratio variables
-      real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
-      real*8  KCOEF(4)
-      integer NVAR(L_LEVELS)
-
-!     temporary variables for multiple aerosol calculation
-      real*8 atemp, btemp
-
-!     variables for k in units m^-1
-      real*8 rho, dz(L_LEVELS)
-
-      integer igas
-
-      !--- Kasting's CIA ----------------------------------------
-      !real*8, parameter :: Ci(L_NSPECTI)=[                         &
-      !     3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7,    &
-      !     5.4E-7, 1.6E-6, 0.0,                               &
-      !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,            & 
-      !     0.0, 4.0E-7, 4.0E-6, 1.4E-5,    &
-      !     1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ] 
-      !real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9,             &
-      !     -1.7, -1.7, -1.7, -1.7, -1.7, -1.7,                &
-      !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-      !     -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
-      !----------------------------------------------------------
-
-!=======================================================================
-!     Determine the total gas opacity throughout the column, for each
-!     spectral interval, NW, and each Gauss point, NG.
-
-      taugsurf(:,:) = 0.0
-      dpr(:)        = 0.0
-      lkcoef(:,:)   = 0.0
-
-      do K=2,L_LEVELS
-         DPR(k) = PLEV(K)-PLEV(K-1)
-
-         !--- Kasting's CIA ----------------------------------------
-         !dz(k)=dpr(k)*189.02*TMID(K)/(0.03720*PMID(K))
-         ! this is CO2 path length (in cm) as written by Francois
-         ! delta_z = delta_p * R_specific * T / (g * P)
-         ! But Kasting states that W is in units of _atmosphere_ cm
-         ! So we do
-         !dz(k)=dz(k)*(PMID(K)/1013.25)
-         !dz(k)=dz(k)/100.0 ! in m for SI calc
-         !----------------------------------------------------------
-
-         ! if we have continuum opacities, we need dz
-         if(kastprof)then
-            dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
-            U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
-         else
-            dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
-            U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
-         endif
-
-         call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
-              LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
-         
-         do LK=1,4
-            LKCOEF(K,LK) = LCOEF(LK)
-         end do
-
-
-         DO NW=1,L_NSPECTI
-            do iaer=1,naerkind
-               TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
-            end do
-         END DO
-      end do                    ! levels
-
-      do K=2,L_LEVELS
-         do nw=1,L_NSPECTI
- 
-            DCONT = 0.0 ! continuum absorption
-
-	    if(Continuum.and.(.not.graybody))then
-            ! include continua if necessary
-               wn_cont = dble(wnoi(nw))
-               T_cont  = dble(TMID(k))
-               do igas=1,ngasmx
-
-                  if(gfrac(igas).eq.-1)then ! variable
-                     p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
-                  else
-                     p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
-                  endif
-
-                  dtemp=0.0
-                  if(gnom(igas).eq.'H2_')then
-                     call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
-                  elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
-                     p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
-                     call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-                  endif
-
-                  DCONT = DCONT + dtemp
-               enddo
-
-               DCONT = DCONT*dz(k) 
-            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
-            !----------------------------------------------------------
-
-
-            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) = GASI(MT(K),MP(K),1,NW,NG)
-                  KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
-                  KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
-                  KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
-               else
-
-                  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),NW,NG))
-
-                  KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*   &
-                       (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                 &
-                       GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
-
-                  KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
-                       (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -               &
-                       GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
-
-                  KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*   &
-                       (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                 &
-                       GASI(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) +            &
-                    LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
-               
-               TAUGAS  = U(k)*ANS
-
-               TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
-               !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
-               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!
-
-            end do
-
-!     Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
-!     which holds continuum opacity only
-
-            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!
-
-         end do
-      end do
-
-
-!=======================================================================
-!     Now the full treatment for the layers, where besides the opacity
-!     we need to calculate the scattering albedo and asymmetry factors
-
-      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
-
-      DO NW=1,L_NSPECTI
-         NG = L_NGAUSS
-         DO L=1,L_NLAYRAD
-
-            K              = 2*L+1
-            DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
-
-            atemp = 0.
-            btemp = 0.
-            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)
-            else
-               WBARI(L,nw,ng) = 0.0D0
-               DTAUI(L,NW,NG) = 1.0E-9
-            endif
-
-            if(btemp .GT. 0.0) then
-               cosbi(L,NW,NG) = atemp/btemp
-            else
-               cosbi(L,NW,NG) = 0.0D0
-            end if
-
-         END DO ! L vertical loop
-
-!     Now the other Gauss points, if needed.
-
-         DO NG=1,L_NGAUSS-1
-            IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
-
-               DO L=1,L_NLAYRAD
-                  K              = 2*L+1
-                  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)
-
-                  else
-                     WBARI(L,nw,ng) = 0.0D0
-                     DTAUI(L,NW,NG) = 1.0E-9
-                  endif
-
-                  cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
-               END DO ! L vertical loop
-            END IF
-            
-         END DO                 ! NG Gauss loop
-      END DO                    ! NW spectral loop
-
-!     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
-
-            TAUCUMI(1,NW,NG)=0.0D0
-            DO K=2,L_LEVELS
-               TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
-            END DO
-         END DO                 ! end full gauss loop
-      END DO
-
-      return
-
-
-    end subroutine optci
-
-
-
+  real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+  real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
+  real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
+  real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
+  real*8 PLEV(L_LEVELS)
+  real*8 TLEV(L_LEVELS)
+  real*8 TMID(L_LEVELS), PMID(L_LEVELS)
+  real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+  real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+
+  ! for aerosols
+  real*8  QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
+  real*8  QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
+  real*8  GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
+  real*8  TAUAERO(L_LEVELS+1,NAERKIND)
+  real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
+  real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
+
+  integer L, NW, NG, K, LK, IAER
+  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
+  real*8  ANS, TAUGAS
+  real*8  DPR(L_LEVELS), U(L_LEVELS)
+  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
+
+  real*8 taugsurf(L_NSPECTI,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)
+  integer NVAR(L_LEVELS)
+
+  ! temporary variables for multiple aerosol calculation
+  real*8 atemp, btemp
+
+  ! variables for k in units m^-1
+  real*8 rho, dz(L_LEVELS)
+
+  integer igas, jgas
+
+  !--- Kasting's CIA ----------------------------------------
+  !real*8, parameter :: Ci(L_NSPECTI)=[                         &
+  !     3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7,    &
+  !     5.4E-7, 1.6E-6, 0.0,                               &
+  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,            & 
+  !     0.0, 4.0E-7, 4.0E-6, 1.4E-5,    &
+  !     1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ] 
+  !real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9,             &
+  !     -1.7, -1.7, -1.7, -1.7, -1.7, -1.7,                &
+  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+  !     -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
+  !----------------------------------------------------------
+
+  !=======================================================================
+  !     Determine the total gas opacity throughout the column, for each
+  !     spectral interval, NW, and each Gauss point, NG.
+
+  taugsurf(:,:) = 0.0
+  dpr(:)        = 0.0
+  lkcoef(:,:)   = 0.0
+
+  do K=2,L_LEVELS
+     DPR(k) = PLEV(K)-PLEV(K-1)
+
+     !--- Kasting's CIA ----------------------------------------
+     !dz(k)=dpr(k)*189.02*TMID(K)/(0.03720*PMID(K))
+     ! this is CO2 path length (in cm) as written by Francois
+     ! delta_z = delta_p * R_specific * T / (g * P)
+     ! But Kasting states that W is in units of _atmosphere_ cm
+     ! So we do
+     !dz(k)=dz(k)*(PMID(K)/1013.25)
+     !dz(k)=dz(k)/100.0 ! in m for SI calc
+     !----------------------------------------------------------
+
+     ! if we have continuum opacities, we need dz
+     if(kastprof)then
+        dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K))
+        U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
+     else
+        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
+        U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
+     endif
+
+     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
+          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
+
+     do LK=1,4
+        LKCOEF(K,LK) = LCOEF(LK)
+     end do
+
+
+     DO NW=1,L_NSPECTI
+        do iaer=1,naerkind
+           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
+        end do
+     END DO
+  end do                    ! levels
+
+  do K=2,L_LEVELS
+     do nw=1,L_NSPECTI
+
+        DCONT = 0.0 ! continuum absorption
+
+        if(Continuum.and.(.not.graybody))then
+           ! include continua if necessary
+           wn_cont = dble(wnoi(nw))
+           T_cont  = dble(TMID(k))
+           do igas=1,ngasmx
+
+              if(gfrac(igas).eq.-1)then ! variable
+                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
+              else
+                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+              endif
+
+              dtemp=0.0
+              if(igas.eq.igas_N2)then
+
+                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
+
+              elseif(igas.eq.igas_H2)then
+
+                 ! first do self-induced absorption
+                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+
+                 ! then cross-interactions with other gases
+                 do jgas=1,ngasmx
+                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    dtempc  = 0.0
+                    if(jgas.eq.igas_N2)then 
+                       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.)
+                    endif
+                    dtemp = dtemp + dtempc
+                 enddo
+
+              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
+
+                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
+                 if(H2Ocont_simple)then
+                    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.)
+                 endif
+
+              endif
+
+              DCONT = DCONT + dtemp
+
+           enddo
+
+           ! Oobleck test
+           !rho = PMID(k)*scalep / (TMID(k)*286.99)
+           !if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
+           !   DCONT = rho * 0.125 * 4.6e-4
+           !elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
+           !   DCONT = 1000*dpr(k) * 1.0 * 4.6e-4 / g
+           !   DCONT = rho * 1.0 * 4.6e-4
+           !elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
+           !   DCONT = rho * 0.125 * 4.6e-4
+           !endif
+
+           DCONT = DCONT*dz(k)
+
+        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
+        !----------------------------------------------------------
+
+        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) = GASI(MT(K),MP(K),1,NW,NG)
+              KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
+              KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
+              KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
+           else
+
+              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),NW,NG))
+
+              KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*   &
+                   (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                 &
+                   GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
+
+              KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
+                   (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -               &
+                   GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
+
+              KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*   &
+                   (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                 &
+                   GASI(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) +            &
+                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
+
+           TAUGAS  = U(k)*ANS
+
+           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!
+
+        end do
+
+        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
+        ! which holds continuum opacity only
+
+        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!
+
+     end do
+  end do
+
+
+  !=======================================================================
+  !     Now the full treatment for the layers, where besides the opacity
+  !     we need to calculate the scattering albedo and asymmetry factors
+
+  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
+
+  DO NW=1,L_NSPECTI
+     NG = L_NGAUSS
+     DO L=1,L_NLAYRAD
+
+        K              = 2*L+1
+        DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
+
+        atemp = 0.
+        btemp = 0.
+        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)
+        else
+           WBARI(L,nw,ng) = 0.0D0
+           DTAUI(L,NW,NG) = 1.0E-9
+        endif
+
+        if(btemp .GT. 0.0) then
+           cosbi(L,NW,NG) = atemp/btemp
+        else
+           cosbi(L,NW,NG) = 0.0D0
+        end if
+
+     END DO ! L vertical loop
+
+     ! Now the other Gauss points, if needed.
+
+     DO NG=1,L_NGAUSS-1
+        IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
+
+           DO L=1,L_NLAYRAD
+              K              = 2*L+1
+              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)
+
+              else
+                 WBARI(L,nw,ng) = 0.0D0
+                 DTAUI(L,NW,NG) = 1.0E-9
+              endif
+
+              cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
+           END DO ! L vertical loop
+        END IF
+
+     END DO                 ! NG Gauss loop
+  END DO                    ! NW spectral loop
+
+  ! 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
+
+        TAUCUMI(1,NW,NG)=0.0D0
+        DO K=2,L_LEVELS
+           TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
+        END DO
+     END DO                 ! end full gauss loop
+  END DO
+
+  ! be aware when comparing with textbook results 
+  ! (e.g. Pierrehumbert p. 218) that 
+  ! taucumi does not take the <cos theta>=0.5 factor into
+  ! account. It is the optical depth for a vertically 
+  ! ascending ray with angle theta = 0.
+
+  !open(127,file='taucum.out')
+  !do nw=1,L_NSPECTI
+  !   write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
+  !enddo
+  !close(127)
+
+  return
+
+
+end subroutine optci
+
+
+
Index: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 716)
@@ -1,300 +1,324 @@
-      SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
-          QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
-          TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
-
-      use radinc_h
-      use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
-      use gases_h
-
-      implicit none
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Calculates shortwave optical constants at each level.
-!     
-!     Authors
-!     -------
-!     Adapted from the NASA Ames code by R. Wordsworth (2009)
-!     
-!==================================================================
-!     
-!     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL  
-!     IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
-!     LAYER: WBAR, DTAU, COSBAR
-!     LEVEL: TAU
-!     
-!     TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
-!     layer L. NW is spectral wavelength interval, ng the Gauss point index.
-!     
-!     TLEV(L) - Temperature at the layer boundary
-!     PLEV(L) - Pressure at the layer boundary (i.e. level)
-!     GASV(NT,NPS,NW,NG) - Visible k-coefficients 
-!     
-!-------------------------------------------------------------------
+SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
+     QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
+     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
+
+  use radinc_h
+  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
+  use gases_h
+
+  implicit none
+
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Calculates shortwave optical constants at each level.
+  !     
+  !     Authors
+  !     -------
+  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
+  !     
+  !==================================================================
+  !     
+  !     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL  
+  !     IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
+  !     LAYER: WBAR, DTAU, COSBAR
+  !     LEVEL: TAU
+  !     
+  !     TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
+  !     layer L. NW is spectral wavelength interval, ng the Gauss point index.
+  !     
+  !     TLEV(L) - Temperature at the layer boundary
+  !     PLEV(L) - Pressure at the layer boundary (i.e. level)
+  !     GASV(NT,NPS,NW,NG) - Visible k-coefficients 
+  !     
+  !-------------------------------------------------------------------
 
 #include "callkeys.h"
 #include "comcstfi.h"
 
-      real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
-      real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
-      real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
-      real*8 PLEV(L_LEVELS)
-      real*8 TMID(L_LEVELS), PMID(L_LEVELS)
-      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
-      integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
-      real*8  ANS, TAUGAS
-      real*8  TRAY(L_LEVELS,L_NSPECTV)
-      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
-
-!     Variable species mixing ratio variables
-      real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
-      real*8 KCOEF(4)
-      integer NVAR(L_LEVELS)
-
-!     temporary variables for multiple aerosol calculation
-      real*8 atemp, btemp, ctemp
-
-!     variables for k in units m^-1
-      double precision wn_cont, p_cont, p_air, T_cont, dtemp
-      real*8 dz(L_LEVELS), DCONT
-
-      integer igas
-
-!=======================================================================
-!     Determine the total gas opacity throughout the column, for each
-!     spectral interval, NW, and each Gauss point, NG.
-!     Calculate the continuum opacities, i.e., those that do not depend on
-!     NG, the Gauss index.
-
-      taugsurf(:,:) = 0.0
-      dpr(:)        = 0.0
-      lkcoef(:,:)   = 0.0
-
-      do K=2,L_LEVELS
-         DPR(k) = PLEV(K)-PLEV(K-1)
-
-         ! if we have continuum opacities, we need dz
-         if(kastprof)then
-            dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
-            U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
-         else
-            dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
-            U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
-         endif
-
-         call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
-              LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
-
-         do LK=1,4
-            LKCOEF(K,LK) = LCOEF(LK)
-         end do
-
-         DO NW=1,L_NSPECTV
-            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
-
-!     we ignore K=1...
-      do K=2,L_LEVELS
-         do NW=1,L_NSPECTV
-
-            TRAYAER = TRAY(K,NW)
-            do iaer=1,naerkind
-               TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
-            end do
-
-            DCONT = 0.0 ! continuum absorption
-
-            if(callgasvis.and.Continuum.and.(.not.graybody))then
-            ! include continua if necessary
-               wn_cont = dble(wnov(nw))
-               T_cont  = dble(TMID(k))
-               do igas=1,ngasmx
-   
-                  if(gfrac(igas).eq.-1)then ! variable
-                     p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
-                  else
-                     p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
-                  endif
-
-                  dtemp=0.0
-                  if(gnom(igas).eq.'H2_')then 
-                     call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
-                  elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
-                     p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
-                     call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-                  endif
-
-                  DCONT = DCONT + dtemp
-
-               enddo
-
-               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
-                  KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
-                  KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
-                  KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
-                  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) -                  &
-                       GASV(MT(K),MP(K),NVAR(K),NW,NG))
-
-                  KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*  &
-                       (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                &
-                       GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
-
-                  KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
-                       (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -              &
-                       GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
-
-                  KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
-                       (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) +           &
-                    LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
-
-               TAUGAS          = U(k)*ANS
-
-
-               !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
-               TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
-               DTAUKV(K,nw,ng) = TAUGAS + TRAYAER  & ! TRAYAER includes all scattering contributions
-                                 + DCONT             ! for continuum absorption
-
-            end do
-
-
-!     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
-            end do ! a bug was here!
-
-         end do
-      end do
-
-
-!=======================================================================
-!     Now the full treatment for the layers, where besides the opacity
-!     we need to calculate the scattering albedo and asymmetry factors
-
-      DO NW=1,L_NSPECTV
-         DO K=2,L_LEVELS
-            do iaer=1,naerkind
-              TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
-            end do
-         ENDDO
-      ENDDO
-      
-
-      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 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
-
-               COSBV(L,NW,NG) = atemp/btemp
-               WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
-
-            END DO
-
-!     No vertical averaging on bottom layer
-
-            L              = L_NLAYRAD
-            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 point loop
-      END DO                    ! NW spectral loop
-
-
-
-!     Total extinction optical depths
-
-      DO NW=1,L_NSPECTV        
-         DO NG=1,L_NGAUSS       ! full gauss loop
-            TAUV(1,NW,NG)=0.0D0
-            DO L=1,L_NLAYRAD
-               TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
-            END DO
-
-            TAUCUMV(1,NW,NG)=0.0D0
-            DO K=2,L_LEVELS
-               TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
-            END DO
-         END DO                 ! end full gauss loop
-      END DO
-
-
-      RETURN
-    END SUBROUTINE OPTCV
+  real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+  real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
+  real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
+  real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
+  real*8 PLEV(L_LEVELS)
+  real*8 TMID(L_LEVELS), PMID(L_LEVELS)
+  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
+  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
+  real*8  ANS, TAUGAS
+  real*8  TRAY(L_LEVELS,L_NSPECTV)
+  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
+
+  ! variable species mixing ratio variables
+  real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
+  real*8 KCOEF(4)
+  integer NVAR(L_LEVELS)
+
+  ! temporary variables for multiple aerosol calculation
+  real*8 atemp, btemp, ctemp
+
+  ! 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
+
+  integer igas, jgas
+
+  !=======================================================================
+  !     Determine the total gas opacity throughout the column, for each
+  !     spectral interval, NW, and each Gauss point, NG.
+  !     Calculate the continuum opacities, i.e., those that do not depend on
+  !     NG, the Gauss index.
+
+  taugsurf(:,:) = 0.0
+  dpr(:)        = 0.0
+  lkcoef(:,:)   = 0.0
+
+  do K=2,L_LEVELS
+     DPR(k) = PLEV(K)-PLEV(K-1)
+
+     ! if we have continuum opacities, we need dz
+     if(kastprof)then
+        dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
+        U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
+     else
+        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
+        U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
+     endif
+
+     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
+          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
+
+     do LK=1,4
+        LKCOEF(K,LK) = LCOEF(LK)
+     end do
+
+     DO NW=1,L_NSPECTV
+        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
+
+  !     we ignore K=1...
+  do K=2,L_LEVELS
+     do NW=1,L_NSPECTV
+
+        TRAYAER = TRAY(K,NW)
+        do iaer=1,naerkind
+           TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
+        end do
+
+        DCONT = 0.0 ! continuum absorption
+
+        if(callgasvis.and.continuum.and.(.not.graybody))then
+           ! include continua if necessary
+           wn_cont = dble(wnov(nw))
+           T_cont  = dble(TMID(k))
+           do igas=1,ngasmx
+
+              if(gfrac(igas).eq.-1)then ! variable
+                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
+              else
+                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+              endif
+
+              dtemp=0.0
+              if(igas.eq.igas_N2)then
+
+                 !call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
+
+              elseif(igas.eq.igas_H2)then
+
+                 ! first do self-induced absorption
+                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+
+                 ! 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.)
+                       ! should be irrelevant in the visible
+                    elseif(jgas.eq.igas_He)then 
+                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
+                    endif
+                 enddo
+
+              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
+
+                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
+                 if(H2Ocont_simple)then
+                    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.)
+                 endif
+
+              endif
+
+              DCONT = DCONT + dtemp
+
+           enddo
+
+           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
+              KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
+              KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
+              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
+              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) -                  &
+                   GASV(MT(K),MP(K),NVAR(K),NW,NG))
+
+              KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*  &
+                   (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                &
+                   GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
+
+              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
+                   (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -              &
+                   GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
+
+              KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
+                   (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) +           &
+                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
+
+           TAUGAS = U(k)*ANS
+
+           !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
+           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
+           DTAUKV(K,nw,ng) = TAUGAS + TRAYAER  & ! TRAYAER includes all scattering contributions
+                + DCONT             ! for continuum absorption
+
+        end do
+
+
+        !     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
+        end do ! a bug was here!
+
+     end do
+  end do
+
+
+  !=======================================================================
+  !     Now the full treatment for the layers, where besides the opacity
+  !     we need to calculate the scattering albedo and asymmetry factors
+
+  DO NW=1,L_NSPECTV
+     DO K=2,L_LEVELS
+        do iaer=1,naerkind
+           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
+        end do
+     ENDDO
+  ENDDO
+
+
+  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 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
+
+           COSBV(L,NW,NG) = atemp/btemp
+           WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
+
+        END DO
+
+        !     No vertical averaging on bottom layer
+
+        L              = L_NLAYRAD
+        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 point loop
+  END DO                    ! NW spectral loop
+
+
+
+  ! Total extinction optical depths
+
+  DO NW=1,L_NSPECTV        
+     DO NG=1,L_NGAUSS       ! full gauss loop
+        TAUV(1,NW,NG)=0.0D0
+        DO L=1,L_NLAYRAD
+           TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
+        END DO
+
+        TAUCUMV(1,NW,NG)=0.0D0
+        DO K=2,L_LEVELS
+           TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
+        END DO
+     END DO                 ! end full gauss loop
+  END DO
+
+
+  RETURN
+END SUBROUTINE OPTCV
Index: trunk/LMDZ.GENERIC/libf/phystd/params_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/params_h.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/params_h.F90	(revision 716)
@@ -6,6 +6,6 @@
       implicit none
 
-      double precision rc,m_n,m_v,rmn
-      save m_n,m_v,rmn
+      double precision rc,cp_n,m_n,m_v,rmn
+      save cp_n,m_n,m_v,rmn
 
       logical ideal_v
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 716)
@@ -100,4 +100,5 @@
 !           To F90: R. Wordsworth (2010)
 !           New turbulent diffusion scheme: J. Leconte (2012)
+!           Loops converted to F90 matrix format: J. Leconte (2012)
 !     
 !==================================================================
@@ -690,6 +691,6 @@
               if(CLFvarying)then
 
-                 !!! ---> PROBLEMS WITH ALLOCATED ARRAYS 
-                 !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
+                 ! ---> PROBLEMS WITH ALLOCATED ARRAYS 
+                 ! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
                  clearsky=.true.
                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
@@ -700,5 +701,5 @@
                       reffrad,tau_col1,cloudfrac,totcloudfrac,             &
                       clearsky,firstcall,lastcall)
-                 clearsky = .false.  !! just in case.     
+                 clearsky = .false.  ! just in case.     
 
                  ! Sum the fluxes and heating rates from cloudy/clear cases
@@ -1866,8 +1867,8 @@
       icount=icount+1
 
-      !!! DEALLOCATE STUFF
+      ! deallocate gas variables
       if (lastcall) then
-        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )    !! this was allocated in su_gases.F90
-        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  !! this was allocated in su_gases.F90
+        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
+        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
       endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 716)
@@ -90,8 +90,8 @@
       integer, parameter :: nsizemax = 60
 
-      character (len=20) :: corrkdir
+      character (len=100) :: corrkdir
       save corrkdir
 
-      character (len=30) :: banddir
+      character (len=100) :: banddir
       save banddir
 
Index: trunk/LMDZ.GENERIC/libf/phystd/setspi.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/setspi.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/setspi.F90	(revision 716)
@@ -36,6 +36,6 @@
 
       character(len=30)  :: temp1
-      character(len=100) :: file_id
-      character(len=100) :: file_path
+      character(len=200) :: file_id
+      character(len=200) :: file_path
 
 !     C1 and C2 values from Goody and Yung (2nd edition)  MKS units
Index: trunk/LMDZ.GENERIC/libf/phystd/setspv.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/setspv.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/setspv.F90	(revision 716)
@@ -38,6 +38,6 @@
 
       character(len=30)  :: temp1
-      character(len=100) :: file_id
-      character(len=100) :: file_path
+      character(len=200) :: file_id
+      character(len=200) :: file_path
 
       real*8 :: lastband(2)
Index: trunk/LMDZ.GENERIC/libf/phystd/su_gases.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/su_gases.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/su_gases.F90	(revision 716)
@@ -5,68 +5,115 @@
   implicit none
 
-       integer igas, ierr
+  integer igas, ierr, count
 
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Load atmospheric composition info
-!     
-!     Authors
-!     -------
-!     R. Wordsworth (2011)
-!     
-!==================================================================
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Load atmospheric composition info
+  !     
+  !     Authors
+  !     -------
+  !     R. Wordsworth (2011)
+  !     Allocatable arrays by A. Spiga (2011)
+  !     
+  !==================================================================
 
-       ! load gas names from file 'gases.def'
-       open(90,file='gases.def',status='old',form='formatted',iostat=ierr)
-       if (ierr.eq.0) then
-          write(*,*) "sugases.F90: reading file gases.def"
-          read(90,*)
-          read(90,*,iostat=ierr) ngasmx
-          if (ierr.ne.0) then
-             write(*,*) "sugases.F90: error reading number of gases"
-             write(*,*) "   (first line of gases.def) "
-             call abort
-          endif
+  ! load gas names from file 'gases.def'
+  open(90,file='gases.def',status='old',form='formatted',iostat=ierr)
+  if (ierr.eq.0) then
+     write(*,*) "sugases.F90: reading file gases.def"
+     read(90,*)
+     read(90,*,iostat=ierr) ngasmx
+     if (ierr.ne.0) then
+        write(*,*) "sugases.F90: error reading number of gases"
+        write(*,*) "   (first line of gases.def) "
+        call abort
+     endif
 
-          PRINT *, "OK I HAVE ",ngasmx, " GASES. NOW I ALLOCATE ARRAYS AND READ NAMES AND FRAC."
+     print*,ngasmx, " gases found in gases.def. Allocating names and molar fractions..."
 
-          IF ( .NOT. ALLOCATED( gnom ) ) ALLOCATE( gnom( ngasmx ) )
-          do igas=1,ngasmx
-             read(90,*,iostat=ierr) gnom(igas)
-             if (ierr.ne.0) then
-                write(*,*) 'sugases.F90: error reading gas names in gases.def...'
-                call abort
-             endif
-          enddo                  !of do igas=1,ngasmx
-          
-          vgas=0
-          IF ( .NOT. ALLOCATED( gfrac ) ) ALLOCATE( gfrac( ngasmx ) )
-          do igas=1,ngasmx
-             read(90,*,iostat=ierr) gfrac(igas)
-             if (ierr.ne.0) then
-                write(*,*) 'sugases.F90: error reading gas molar fractions in gases.def...'
-                call abort
-             endif
+     if (.not.allocated(gnom)) allocate(gnom(ngasmx))
+     do igas=1,ngasmx
+        read(90,*,iostat=ierr) gnom(igas)
+        if (ierr.ne.0) then
+           write(*,*) 'sugases.F90: error reading gas names in gases.def...'
+           call abort
+        endif
+     enddo                  !of do igas=1,ngasmx
 
-             ! find variable gas (if any)
-             if(gfrac(igas).eq.-1.0)then 
-                if(vgas.eq.0)then
-                   vgas=igas
-                else
-                   print*,'You seem to be choosing two variable gases'
-                   print*,'Check that gases.def is correct'
-                   call abort
-                endif
-             endif
-             
-          enddo                  !of do igas=1,ngasmx
+     vgas=0
+     if(.not.allocated(gfrac)) allocate(gfrac(ngasmx))
+     do igas=1,ngasmx
+        read(90,*,iostat=ierr) gfrac(igas)
+        if (ierr.ne.0) then
+           write(*,*) 'sugases.F90: error reading gas molar fractions in gases.def...'
+           call abort
+        endif
 
-       else
-          write(*,*) 'Cannot find required file "gases.def"'
-          call abort
-       endif
-       close(90)
+        ! find variable gas (if any)
+        if(gfrac(igas).eq.-1.0)then 
+           if(vgas.eq.0)then
+              vgas=igas
+           else
+              print*,'You seem to be choosing two variable gases'
+              print*,'Check that gases.def is correct'
+              call abort
+           endif
+        endif
+
+     enddo                  !of do igas=1,ngasmx
+
+
+     ! assign the 'igas_X' labels
+     count=0
+     do igas=1,ngasmx
+        if (gnom(igas).eq."H2_") then
+           igas_H2=igas
+           count=count+1
+        elseif (gnom(igas).eq."He_") then
+           igas_He=igas
+           count=count+1
+        elseif (gnom(igas).eq."H2O") then
+           igas_H2O=igas
+           count=count+1
+        elseif (gnom(igas).eq."CO2") then
+           igas_CO2=igas
+           count=count+1
+        elseif (gnom(igas).eq."CO_") then
+           igas_CO=igas
+           count=count+1
+        elseif (gnom(igas).eq."N2_") then
+           igas_N2=igas
+           count=count+1
+        elseif (gnom(igas).eq."O2_") then
+           igas_O2=igas
+           count=count+1
+        elseif (gnom(igas).eq."SO2") then
+           igas_SO2=igas
+           count=count+1
+        elseif (gnom(igas).eq."H2S") then
+           igas_H2S=igas
+           count=count+1
+        elseif (gnom(igas).eq."CH4") then
+           igas_CH4=igas
+           count=count+1
+        elseif (gnom(igas).eq."NH3") then
+           igas_NH3=igas
+           count=count+1
+        endif
+     enddo
+
+     if(count.ne.ngasmx)then
+        print*,'Mismatch between ngas and number of recognised gases in sugas_corrk.F90.'
+        print*,'Either we haven`t managed to assign all the gases, or there are duplicates.'
+        print*,'Please try again.'
+     endif
+
+  else
+     write(*,*) 'Cannot find required file "gases.def"'
+     call abort
+  endif
+  close(90)
 
 end subroutine su_gases
Index: trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 715)
+++ trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 716)
@@ -13,4 +13,5 @@
 !     Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009)
 !     Added double gray case by Jeremy Leconte (2012)
+!     New HITRAN continuum data section by RW (2012)
 !
 !     Summary
@@ -25,4 +26,5 @@
       use radcommon_h, only : wrefvar
       use datafile_mod, only: datadir
+
       use gases_h
       use ioipsl_getincom 
@@ -31,4 +33,5 @@
 #include "callkeys.h"
 #include "comcstfi.h"
+
 !==================================================================
 
@@ -38,8 +41,8 @@
       integer L_NGAUSScheck
 
-      character(len=100) :: file_id
-      character(len=300) :: file_path
-
-      !! ALLOCATABLE ARRAYS -- AS 12/2011
+      character(len=200) :: file_id
+      character(len=500) :: file_path
+
+      ! ALLOCATABLE ARRAYS -- AS 12/2011
       REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi8, gasv8
       character*3,allocatable,DIMENSION(:) :: gastype ! for check with gnom
@@ -48,5 +51,5 @@
       real kappa_IR, kappa_VI
 
-      integer ngas, igas
+      integer ngas, igas, jgas
 
       double precision testcont ! for continuum absorption initialisation
@@ -84,10 +87,5 @@
                 '] has no variable species, exiting.'
          call abort
-      elseif(ngas.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
-         print*,'You have varactive and varfixed=.false. and the database [', &
-                     corrkdir(1:LEN_TRIM(corrkdir)),           &
-                '] has a variable species.'
-         call abort
-      elseif(ngas.gt.4 .or. ngas.lt.1)then
+      elseif(ngas.gt.5 .or. ngas.lt.1)then
          print*,ngas,' species in database [',               &
                      corrkdir(1:LEN_TRIM(corrkdir)),           &
@@ -96,15 +94,13 @@
       endif 
 
-      if(ngas.gt.3)then
-         print*,'ngas>3, EXPERIMENTAL!'
-      endif
-
+      ! dynamically allocate gastype and read from Q.dat
       IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) )
-      do n=1,ngas
-         read(111,*) gastype(n)
-         print*,'Gas ',n,' is ',gastype(n)
+
+      do igas=1,ngas
+         read(111,*) gastype(igas)
+         print*,'Gas ',igas,' is ',gastype(igas)
       enddo
 
-      ! GET array size, load the coefficients
+      ! get array size, load the coefficients
       open(111,file=TRIM(file_path),form='formatted')
       read(111,*) L_REFVAR
@@ -113,10 +109,17 @@
       close(111)
 
+      if(L_REFVAR.gt.1 .and. (.not.varactive) .and. (.not.varfixed))then
+         print*,'You have varactive and varfixed=.false. and the database [', &
+                     corrkdir(1:LEN_TRIM(corrkdir)),           &
+                '] has a variable species.'
+         call abort
+      endif
+
       ! Check that gastype and gnom match
-      do n=1,ngas
-         print*,'Gas ',n,' is ',gnom(n)
-         if(gnom(n).ne.gastype(n))then
-            print*,'Name of a gas in radiative transfer data (',gastype(n),') does not ', &
-                 'match that in gases.def (',gnom(n),'), exiting. You should compare ', &
+      do igas=1,ngas
+         print*,'Gas ',igas,' is ',gnom(igas)
+         if(gnom(igas).ne.gastype(igas))then
+            print*,'Name of a gas in radiative transfer data (',gastype(igas),') does not ', &
+                 'match that in gases.def (',gnom(igas),'), exiting. You should compare ', &
                  'gases.def with Q.dat in your radiative transfer directory.' 
             call abort
@@ -126,5 +129,5 @@
 
       ! display the values
-      print*,'Variable gas mixing ratios:'
+      print*,'Variable gas volume mixing ratios:'
       do n=1,L_REFVAR
          !print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
@@ -191,5 +194,5 @@
       endif
      
-      ! GET array size, load the coefficients
+      ! get array size, load the coefficients
       open(111,file=TRIM(file_path),form='formatted')
       read(111,*) L_NPREF
@@ -200,5 +203,4 @@
       IF( .NOT. ALLOCATED( pfgasref ) ) ALLOCATE( PFGASREF(L_PINT) )
 
-
       ! display the values
       print*,'Correlated-k pressure grid (mBar):'
@@ -219,5 +221,5 @@
       end do
       pfgasref(L_PINT) = pgasref(L_NPREF)
-      print*,'Warning, pfgasref needs generalisation to uneven grids!!'
+      ! Warning, pfgasref needs generalisation to uneven grids!
 
 !-----------------------------------------------------------------------
@@ -239,5 +241,5 @@
       endif
 
-      ! GET array size, load the coefficients
+      ! get array size, load the coefficients
       open(111,file=TRIM(file_path),form='formatted')
       read(111,*) L_NTREF
@@ -256,15 +258,16 @@
       tgasmax = tgasref(L_NTREF)
 
-
 !-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-! ALLOCATE THE MULTIDIM ARRAYS IN radcommon_h
-      PRINT *, L_NTREF,L_NPREF,L_REFVAR 
+! allocate the multidimensional arrays in radcommon_h
+
       IF( .NOT. ALLOCATED( gasi8 ) ) ALLOCATE( gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS) )
       IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
       IF( .NOT. ALLOCATED( gasi ) ) ALLOCATE( gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS) )
       IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
-!-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
+
+      ! display the values
+      print*,''
+      print*,'Correlated-k matrix size:' 
+      print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']' 
 
 !=======================================================================
@@ -288,16 +291,18 @@
          read(111,*) gasv8
          close(111)
-	 
+
       else if (callgasvis.and.graybody) then
-!    constant absorption coefficient in visible
+         ! constant absorption coefficient in visible
          write(*,*)"constant absorption coefficient in visible:"
-         kappa_VI=-100000.
+         kappa_VI = -100000.
          call getin("kappa_VI",kappa_VI)
          write(*,*)" kappa_VI = ",kappa_VI
-	 kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule         
-         gasv8(:,:,:,:,:)=kappa_VI
+	 kappa_VI = kappa_VI * 1.e4 * mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule         
+         gasv8(:,:,:,:,:) = kappa_VI
+
       else
          print*,'Visible corrk gaseous absorption is set to zero.'
-         gasv8(:,:,:,:,:)=0.0
+         gasv8(:,:,:,:,:) = 0.0
+
       endif
 
@@ -319,25 +324,57 @@
          read(111,*) gasi8
          close(111)
+     
+         ! 'fzero' is a currently unused feature that allows optimisation
+         ! of the radiative transfer by neglecting bands where absorption
+         ! is close to zero. As it could be useful in the future, this 
+         ! section of the code has been kept commented and not erased.
+         ! RW 7/3/12.
+
+         do nw=1,L_NSPECTI
+            fzeroi(nw) = 0.
+!            do nt=1,L_NTREF
+!               do np=1,L_NPREF
+!                  do nh=1,L_REFVAR
+!                     do ng = 1,L_NGAUSS
+!                        if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
+!                           fzeroi(nw)=fzeroi(nw)+1.
+!                        endif
+!                     end do
+!                  end do
+!               end do
+!            end do
+!            fzeroi(nw)=fzeroi(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
+         end do
+
+         do nw=1,L_NSPECTV
+            fzerov(nw) = 0.
+!            do nt=1,L_NTREF
+!               do np=1,L_NPREF
+!                  do nh=1,L_REFVAR
+!                     do ng = 1,L_NGAUSS
+!                        if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
+!                           fzerov(nw)=fzerov(nw)+1.
+!                        endif
+!                     end do
+!                  end do
+!               end do
+!            end do
+!            fzerov(nw)=fzerov(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
+         end do
 
       else if (graybody) then
-!     constant absorption coefficient in IR
+         ! constant absorption coefficient in IR
          write(*,*)"constant absorption coefficient in InfraRed:"
-         kappa_IR=-100000.
+         kappa_IR = -100000.
          call getin("kappa_IR",kappa_IR)
-         write(*,*)" kappa_IR = ",kappa_IR	 
-         kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule       
-         gasi8(:,:,:,:,:)=kappa_IR
+         write(*,*)" kappa_IR = ",kappa_IR
+         kappa_IR = kappa_IR * 1.e4 * mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule       
+         gasi8(:,:,:,:,:) = kappa_IR
+
       else
          print*,'Infrared corrk gaseous absorption is set to zero.'
-         gasi8(:,:,:,:,:)=0.0
-      endif
-
-      do nw=1,L_NSPECTI
-         fzeroi(nw) = 0.
-      end do
-
-      do nw=1,L_NSPECTV
-         fzerov(nw) = 0.
-      end do
+         gasi8(:,:,:,:,:) = 0.0
+
+      endif
 
 
@@ -512,14 +549,46 @@
 
 
+!=======================================================================
+!     Initialise the continuum absorption data
 
       do igas=1,ngasmx
-         if(gnom(igas).eq.'H2_')then 
+
+         if(gnom(igas).eq.'N2_')then 
+
+            call interpolateN2N2(100.D+0,250.D+0,17500.D+0,testcont,.true.)
+
+         elseif(gnom(igas).eq.'H2_')then
+
+            ! first do self-induced absorption
             call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.)
-         elseif(gnom(igas).eq.'H2O')then 
-            call interpolateH2Ocont(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
+            ! then cross-interactions with other gases
+            do jgas=1,ngasmx
+               if(gnom(jgas).eq.'N2_')then 
+                  call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.)
+               elseif(gnom(jgas).eq.'He_')then 
+                  call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.)
+               endif
+            enddo
+
+         elseif(gnom(igas).eq.'H2O')then
+
+            ! H2O is special 
+            if(H2Ocont_simple)then
+               call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
+            else
+               call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
+            endif
+
          endif
+
       enddo
 
-      !!! DEALLOCATE LOCAL ARRAYS
+      print*,'----------------------------------------------------'
+      print*,'And that`s all we have. It`s possible that other'
+      print*,'continuum absorption may be present, but if it is we'
+      print*,'don`t yet have data for it...'
+      print*,''
+
+!     Deallocate local arrays
       IF( ALLOCATED( gasi8 ) ) DEALLOCATE( gasi8 )
       IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
