Index: trunk/LMDZ.GENERIC/utilities/co2_cont.for
===================================================================
--- trunk/LMDZ.GENERIC/utilities/co2_cont.for	(revision 1507)
+++ 	(revision )
@@ -1,346 +1,0 @@
-c     This program computes RT CIA spectra of CO2 at low density limit 
-C     at temperatures between 200 and 800K.
-C     
-C     ============================================================================
-C     Copyright (C) 1997   Marcin Gruszka 
-C     ============================================================================
-C     
-C     Copyright Notice:
-C     You may use this program for your scientific applications,
-C     but please do not	distribute it yourself.
-C     Direct all your inquires to the  author: e-mail aborysow@stella.nbi.dk
-C     Final request: If you publish  your work which benefited from this
-C     program, please acknowledge using this program and quote 
-C     the original paper describing the procedure:
-C     M. Gruszka and A. Borysow,
-C     Roto-translational collision-induced absorption of CO2
-C     for the Atmosphere of Venus at Frequencies from 0 to 250 cm$^{-1}$
-C     and at temperature from 200K to 800K,
-C     Icarus, vol. 129, pp. 172-177, (1997). 
-
-C     
-C     The paper (Icarus) describing this computer model, is based on  original paper by
-C     M. Gruszka and A. Borysow,
-C     "Computer Simulation of the Far Infrared Collision-induced absorption spectra 
-C     of gaseous CO2", Mol. Phys., vol. 93, pp. 1007-1016, 1998. 
-C     ============================================================================
-C     
-C     NOTE: The program gives data at spacing corresponding to
-c     Freq(max)-Freq(min)/(npoint-1),
-c     Since reasons are unknown (author's secret), it is recommended to
-c     request one more point that neccesary and get the right freq. step
-c     ==================================================================
-
-c     
-c     The following input parametres are required (*):
-c     ------------------------------------------------
-c     (*) some input parameters are Fortran REAL numbers
-c     and should be entered with a decimal point,
-c     not like INTEGER!
-c     
-c     1. Temperature: (in K) (REAL)
-c     -the model is restricted to the temperature range from
-c     200 K to 800 K.
-c     
-c     2. Frequency limits: (in cm^-1) (REAL,REAL)
-c     -the model has been designed to reproduce the MD and the
-c     experimental data up to 250 cm^-1. There is no restriction
-c     on requesting a higher value, but the accuracy may decrease 
-c     significantly.
-c     
-c     3. Number of points in the output file: (less than Npmx) (INTEGER)
-c     -the frequency grid is controled by the number of points
-c     in the output file.
-c     (i.e. if the upper limit = 250.0 and the number of points = 50.0,
-c     the spectrum will be computed every 5 cm^-1 )
-c     
-c     4. The output in automaticaly directed to output.dat file.
-c     
-c     STRUCTURE OF THE PROGRAM:
-c     -------------------------------------------------------------- 
-c     The MAIN program consists of three modules:
-c     
-c     1. getdat - provides the user interface, 
-c     in this subroutine the Temperature, Frequency Limits and
-c     the number of points are read from the keyboard.
-c     
-c     2. getspc - the main subroutine, where the spectrum is computed.
-c     The three input parameters (temp, range and n; all REAL numbers)
-c     are transmited by value and the absorption coefficient abcoef(..) 
-c     is returned. The default lenght of abcoef(..) is Npmx and if 
-c     n < Npmx, the empty places are set to 0. If a denser grid is
-c     required please change the n appropriaetly in the data declaration
-c     of the MAIN program and inside the getspc.
-c     
-c     3. putdat - generates the output file.
-c     The frequency (in cm^-1) and 
-c     the absorption coefficient (in cm^-1 amagat^-2) are written
-c     to the 'output.dat' file in the followinf format:
-c     format(f10.3,e20.7)
-c     
-c===========================================================================
-c     
-c     
-c     ---------------------------------------------------------------------------------------------------------------
-c     2008/10/15
-c     
-c     Modifications of the code by V. Eymet as follows:
-c     + only the "getspc" routine remains
-c     + inputs: temp (temperature in K), wn (wavenumber in cm^-1); output: abcoef (absorption coefficient, in cm^-1 amagat^-2)
-c     ---------------------------------------------------------------------------------------------------------------
-c
-c
-c
-c     ---------------------------------------------------------------------------------------------------------------
-c     2009/2/28
-c     
-c     Modifications of the code by R. Wordsworth as follows:
-c     + "baranov" routine added to calculate CIA dimer spectrum between 1100 and 1600 cm^-1
-c     + this uses recorded data rather than an analytical function, so bilinear interpolation is needed
-c     + inputs and outputs are identical to those of getspc
-c     ---------------------------------------------------------------------------------------------------------------
-
-
-
-      subroutine getspc(temp,wn,abcoef) ! computes the spectrum
-      implicit none
-      real*8 temp               !temperature
-      real*8 wn        !wavenumber
-      real*8 abcoef       !absorption coefficient
-c     -------------------------------------------------------
-c     The A, B and C parameters as given in the paper:
-      integer ndeg
-      parameter (ndeg=3)
-      real*8 ah(ndeg),bh(ndeg),at(ndeg),bt(ndeg),gam(ndeg)
-c     tau^{L}_{1}:
-      data ah /0.1586914382D+13,-0.9344296879D+01,0.6943966881D+00/
-c     tau^{L}_{2}:
-      data bh /0.1285676961D-12,0.9420973263D+01,-0.7855988401D+00/
-c     tau^{H}_{1}:
-      data at /0.3312598766D-09,0.7285659464D+01,-0.6732642658D+00/
-c     tau^{H}_{2}:
-      data bt /0.1960966173D+09,-0.6834613750D+01,0.5516825232D+00/
-c     gamma_{1}:
-      data gam /0.1059151675D+17,-0.1048630307D+02, 0.7321430968D+00/
-c     **************************************************************
-c     The S - shape function parameters:
-      real*8 w1,w2
-      data w1,w2 /50.0d0, 100.0d0/
-c     --------------------------------------------------------------
-c     local variables:
-c     ----------------------------------
-      real*8 incrmt
-      integer p1,p2,n
-      real*8 a1,b1,a2,b2,gamma
-      real*8 bc1,bc2,bcbc
-      real*8 scon,mtot,spunit,gm0con
-      real*8 icm, frq, lntemp
-      real*8 xk1,x1,x2
-      external xk1
-      incrmt=250.0D+0
-      n=1
-      p1=idint(dnint(w1/incrmt))
-      p2=idint(dnint(w2/incrmt))
-      lntemp=dlog(temp)
-      a1=1.0d0/(ah(1)*dexp(ah(2)*lntemp+ah(3)*lntemp*lntemp))
-      b1=bh(1)*dexp(bh(2)*lntemp+bh(3)*lntemp*lntemp)
-      a2=1.0d0/(at(1)*dexp(at(2)*lntemp+at(3)*lntemp*lntemp))
-      b2=bt(1)*dexp(bt(2)*lntemp+bt(3)*lntemp*lntemp)
-      gamma=gam(1)*dexp(gam(2)*lntemp+gam(3)*lntemp*lntemp)
-      icm=0.1885d0
-      frq=wn
-      x1=b1*dsqrt(a1*a1+icm*icm*frq*frq)
-      bc1=dexp(a1*b1)*a1*XK1(x1)/(a1*a1+icm*icm*frq*frq)
-      x2=b2*dsqrt(a2*a2+icm*icm*frq*frq)
-      bc2=dexp(a2*b2)*a2*XK1(x2)/(a2*a2+icm*icm*frq*frq)
-      if (p1.ge.1) then
-         bcbc=bc1
-      else if (p2.lt.1) then
-         bcbc=bc2
-      else
-         bcbc=dexp( (1.0d0-dble(1-p1)/dble(p2-p1))*dlog(bc1)+
-     &        dble(1-p1)/dble(p2-p1)*dlog(bc2) )
-      endif
-      spunit=1.296917d55
-      gm0con=1.259009d-6
-      scon=spunit/temp
-      mtot=gm0con*temp*gamma*1.0d-56
-      frq=wn
-      abcoef=scon*mtot*bcbc*frq*frq
-      return
-      end
-c     
-c     
-      FUNCTION XK1(X)
-C     MODIFIED BESSEL FUNCTION K1(X) TIMES X
-C     PRECISION IS BETTER THAN 2.2e-7 EVERYWHERE.
-C     ABRAMOWITZ AND S,TEGUN, P.379; TABLES P.417.
-      implicit double precision (a-h,o-z)
-!      IF(X-2.) 10,10,20
-      IF(X-2..le.0.0)then
-! 10   T=(X/3.75)**2
-      T=(X/3.75)**2
-      FI1=X*((((((.00032411*T+.00301532)*T+.02658733)*T+.15084934)
-     1     *T+.51498869)*T+.87890594)*T+.5)
-      T=(X/2.)**2
-      P=(((((-.00004686*T-.00110404)*T-.01919402)*T-.18156897)*T-
-     1     .67278579)*T+.15443144)*T+1.
-      XK1=X*dLOG(X/2)*FI1+P
-!      RETURN
-      else
-! 20   T=2./X
-      T=2./X
-      P=(((((-.00068245*T+.00325614)*T-.00780353)*T+.01504268)*T-
-     1     .03655620)*T+.23498619)*T+1.25331414
-      X=dMIN1(X,330.d0)
-      XK1=dSQRT(X)*dEXP(-X)*P
-      endif
-      RETURN
-      END
-
-
-
-
-c     ----------------------------------------------------------------------------------------------
-c     Added by RDW for inclusion of Baranov (2004) data
-      subroutine baranov(temp,wn,abcoef)
-      ! computes the dimer spectrum
-
-      implicit none
-      real*8 temp               !temperature
-      real*8 wn                 !wavenumber
-      real*8 abcoef             !absorption coefficient
-
-      integer nS,nT
-!      parameter(nS=1713)
-      parameter(nS=3000)
-      parameter(nT=9)
-
-      real*8 wn_arr(nS)
-      real*8 temp_arr(nT)
-      real*8 dim_arr(nS,nT)
-
-      integer pop
-      logical firstcall
-
-      save wn_arr, temp_arr, dim_arr
-
-      character*100 dt_file
-      integer strlen,ios
-      character*100 label
-      label='subroutine baranov'
-
-      dt_file='/san/home/rdword/kspectrum/line_data/CO2_dimer_data_NEW'
-!      dt_file='./data/CO2_dimer_data'
-      open(33,file=dt_file(1:strlen(dt_file)),
-     &     form='unformatted',
-     &     status='old',iostat=ios)
-      if (ios.ne.0) then        ! file not found
-         write(*,*) 'Error from ',label(1:strlen(label)),' :'
-         write(*,*) 'Data file could not be found:'
-         write(*,*) dt_file(1:strlen(dt_file))
-         stop
-      else
-         read(33) wn_arr
-         read(33) temp_arr
-         read(33) dim_arr
-      endif
-      close(33)
-
-
-!	print*,wn_arr(1)!
-!	print*,wn_arr(3000)
-!	print*,temp_arr
-!	print*,dim_arr(1500,1)
-!	stop
-
-
-      call bilinear(wn_arr,temp_arr,nS,nT,dim_arr,wn,temp,abcoef)
-
- 111  continue
-      return
-      end
-c     ----------------------------------------------------------------------------------------------
-
-c     ----------------------------------------------------------------------------------------------
-c     Added by RDW for inclusion of Baranov (2004) data
-      subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x,y,f)
-
-      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_arr(nX)
-      real*8 y_arr(nY)
-      real*8 f2d_arr(nX,nY)
-      
-      integer strlen
-      character*100 label
-      label='subroutine bilinear'
-
-c     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
-c     print*,'x=',x,x_Arr(2),x_arr(nX-2)
-c     stop('Outside CIA dimer wavenumber range!\n')
-         goto 112
-      else
-         
-c     in the x (wavenumber) direction 1st
-         i=1
- 10      if (i.lt.(nX+1)) then ! what passes for a 'while' loop in this sorry language...
-            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(*,*) 'Error from ',label(1:strlen(label)),' :'
-         write(*,*) 'Outside CIA dimer 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
-!         stop
-      else
-c     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)
-      
-c     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)
-      
-c     then in y-direction
-      f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
-      
- 112  continue
-      return
-      end
-c     ----------------------------------------------------------------------------------------------
-      
Index: trunk/LMDZ.GENERIC/utilities/compile
===================================================================
--- trunk/LMDZ.GENERIC/utilities/compile	(revision 1507)
+++ 	(revision )
@@ -1,11 +1,0 @@
-#!/bin/bash
-#gfortran -c -g -fbounds-check generate_kmatrix.F90 strlen.for
-
-gfortran -c generate_kmatrix.F90 co2_cont.for strlen.for
-gfortran -o generate_kmatrix.exe generate_kmatrix.o co2_cont.o strlen.o
-
-gfortran -c make_composition_gen.for strlen.for
-gfortran -o make_composition_gen.exe make_composition_gen.o strlen.o
-
-rm -f *o
-rm -f *~
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/co2_cont.for
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/co2_cont.for	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/co2_cont.for	(revision 1512)
@@ -0,0 +1,346 @@
+c     This program computes RT CIA spectra of CO2 at low density limit 
+C     at temperatures between 200 and 800K.
+C     
+C     ============================================================================
+C     Copyright (C) 1997   Marcin Gruszka 
+C     ============================================================================
+C     
+C     Copyright Notice:
+C     You may use this program for your scientific applications,
+C     but please do not	distribute it yourself.
+C     Direct all your inquires to the  author: e-mail aborysow@stella.nbi.dk
+C     Final request: If you publish  your work which benefited from this
+C     program, please acknowledge using this program and quote 
+C     the original paper describing the procedure:
+C     M. Gruszka and A. Borysow,
+C     Roto-translational collision-induced absorption of CO2
+C     for the Atmosphere of Venus at Frequencies from 0 to 250 cm$^{-1}$
+C     and at temperature from 200K to 800K,
+C     Icarus, vol. 129, pp. 172-177, (1997). 
+
+C     
+C     The paper (Icarus) describing this computer model, is based on  original paper by
+C     M. Gruszka and A. Borysow,
+C     "Computer Simulation of the Far Infrared Collision-induced absorption spectra 
+C     of gaseous CO2", Mol. Phys., vol. 93, pp. 1007-1016, 1998. 
+C     ============================================================================
+C     
+C     NOTE: The program gives data at spacing corresponding to
+c     Freq(max)-Freq(min)/(npoint-1),
+c     Since reasons are unknown (author's secret), it is recommended to
+c     request one more point that neccesary and get the right freq. step
+c     ==================================================================
+
+c     
+c     The following input parametres are required (*):
+c     ------------------------------------------------
+c     (*) some input parameters are Fortran REAL numbers
+c     and should be entered with a decimal point,
+c     not like INTEGER!
+c     
+c     1. Temperature: (in K) (REAL)
+c     -the model is restricted to the temperature range from
+c     200 K to 800 K.
+c     
+c     2. Frequency limits: (in cm^-1) (REAL,REAL)
+c     -the model has been designed to reproduce the MD and the
+c     experimental data up to 250 cm^-1. There is no restriction
+c     on requesting a higher value, but the accuracy may decrease 
+c     significantly.
+c     
+c     3. Number of points in the output file: (less than Npmx) (INTEGER)
+c     -the frequency grid is controled by the number of points
+c     in the output file.
+c     (i.e. if the upper limit = 250.0 and the number of points = 50.0,
+c     the spectrum will be computed every 5 cm^-1 )
+c     
+c     4. The output in automaticaly directed to output.dat file.
+c     
+c     STRUCTURE OF THE PROGRAM:
+c     -------------------------------------------------------------- 
+c     The MAIN program consists of three modules:
+c     
+c     1. getdat - provides the user interface, 
+c     in this subroutine the Temperature, Frequency Limits and
+c     the number of points are read from the keyboard.
+c     
+c     2. getspc - the main subroutine, where the spectrum is computed.
+c     The three input parameters (temp, range and n; all REAL numbers)
+c     are transmited by value and the absorption coefficient abcoef(..) 
+c     is returned. The default lenght of abcoef(..) is Npmx and if 
+c     n < Npmx, the empty places are set to 0. If a denser grid is
+c     required please change the n appropriaetly in the data declaration
+c     of the MAIN program and inside the getspc.
+c     
+c     3. putdat - generates the output file.
+c     The frequency (in cm^-1) and 
+c     the absorption coefficient (in cm^-1 amagat^-2) are written
+c     to the 'output.dat' file in the followinf format:
+c     format(f10.3,e20.7)
+c     
+c===========================================================================
+c     
+c     
+c     ---------------------------------------------------------------------------------------------------------------
+c     2008/10/15
+c     
+c     Modifications of the code by V. Eymet as follows:
+c     + only the "getspc" routine remains
+c     + inputs: temp (temperature in K), wn (wavenumber in cm^-1); output: abcoef (absorption coefficient, in cm^-1 amagat^-2)
+c     ---------------------------------------------------------------------------------------------------------------
+c
+c
+c
+c     ---------------------------------------------------------------------------------------------------------------
+c     2009/2/28
+c     
+c     Modifications of the code by R. Wordsworth as follows:
+c     + "baranov" routine added to calculate CIA dimer spectrum between 1100 and 1600 cm^-1
+c     + this uses recorded data rather than an analytical function, so bilinear interpolation is needed
+c     + inputs and outputs are identical to those of getspc
+c     ---------------------------------------------------------------------------------------------------------------
+
+
+
+      subroutine getspc(temp,wn,abcoef) ! computes the spectrum
+      implicit none
+      real*8 temp               !temperature
+      real*8 wn        !wavenumber
+      real*8 abcoef       !absorption coefficient
+c     -------------------------------------------------------
+c     The A, B and C parameters as given in the paper:
+      integer ndeg
+      parameter (ndeg=3)
+      real*8 ah(ndeg),bh(ndeg),at(ndeg),bt(ndeg),gam(ndeg)
+c     tau^{L}_{1}:
+      data ah /0.1586914382D+13,-0.9344296879D+01,0.6943966881D+00/
+c     tau^{L}_{2}:
+      data bh /0.1285676961D-12,0.9420973263D+01,-0.7855988401D+00/
+c     tau^{H}_{1}:
+      data at /0.3312598766D-09,0.7285659464D+01,-0.6732642658D+00/
+c     tau^{H}_{2}:
+      data bt /0.1960966173D+09,-0.6834613750D+01,0.5516825232D+00/
+c     gamma_{1}:
+      data gam /0.1059151675D+17,-0.1048630307D+02, 0.7321430968D+00/
+c     **************************************************************
+c     The S - shape function parameters:
+      real*8 w1,w2
+      data w1,w2 /50.0d0, 100.0d0/
+c     --------------------------------------------------------------
+c     local variables:
+c     ----------------------------------
+      real*8 incrmt
+      integer p1,p2,n
+      real*8 a1,b1,a2,b2,gamma
+      real*8 bc1,bc2,bcbc
+      real*8 scon,mtot,spunit,gm0con
+      real*8 icm, frq, lntemp
+      real*8 xk1,x1,x2
+      external xk1
+      incrmt=250.0D+0
+      n=1
+      p1=idint(dnint(w1/incrmt))
+      p2=idint(dnint(w2/incrmt))
+      lntemp=dlog(temp)
+      a1=1.0d0/(ah(1)*dexp(ah(2)*lntemp+ah(3)*lntemp*lntemp))
+      b1=bh(1)*dexp(bh(2)*lntemp+bh(3)*lntemp*lntemp)
+      a2=1.0d0/(at(1)*dexp(at(2)*lntemp+at(3)*lntemp*lntemp))
+      b2=bt(1)*dexp(bt(2)*lntemp+bt(3)*lntemp*lntemp)
+      gamma=gam(1)*dexp(gam(2)*lntemp+gam(3)*lntemp*lntemp)
+      icm=0.1885d0
+      frq=wn
+      x1=b1*dsqrt(a1*a1+icm*icm*frq*frq)
+      bc1=dexp(a1*b1)*a1*XK1(x1)/(a1*a1+icm*icm*frq*frq)
+      x2=b2*dsqrt(a2*a2+icm*icm*frq*frq)
+      bc2=dexp(a2*b2)*a2*XK1(x2)/(a2*a2+icm*icm*frq*frq)
+      if (p1.ge.1) then
+         bcbc=bc1
+      else if (p2.lt.1) then
+         bcbc=bc2
+      else
+         bcbc=dexp( (1.0d0-dble(1-p1)/dble(p2-p1))*dlog(bc1)+
+     &        dble(1-p1)/dble(p2-p1)*dlog(bc2) )
+      endif
+      spunit=1.296917d55
+      gm0con=1.259009d-6
+      scon=spunit/temp
+      mtot=gm0con*temp*gamma*1.0d-56
+      frq=wn
+      abcoef=scon*mtot*bcbc*frq*frq
+      return
+      end
+c     
+c     
+      FUNCTION XK1(X)
+C     MODIFIED BESSEL FUNCTION K1(X) TIMES X
+C     PRECISION IS BETTER THAN 2.2e-7 EVERYWHERE.
+C     ABRAMOWITZ AND S,TEGUN, P.379; TABLES P.417.
+      implicit double precision (a-h,o-z)
+!      IF(X-2.) 10,10,20
+      IF(X-2..le.0.0)then
+! 10   T=(X/3.75)**2
+      T=(X/3.75)**2
+      FI1=X*((((((.00032411*T+.00301532)*T+.02658733)*T+.15084934)
+     1     *T+.51498869)*T+.87890594)*T+.5)
+      T=(X/2.)**2
+      P=(((((-.00004686*T-.00110404)*T-.01919402)*T-.18156897)*T-
+     1     .67278579)*T+.15443144)*T+1.
+      XK1=X*dLOG(X/2)*FI1+P
+!      RETURN
+      else
+! 20   T=2./X
+      T=2./X
+      P=(((((-.00068245*T+.00325614)*T-.00780353)*T+.01504268)*T-
+     1     .03655620)*T+.23498619)*T+1.25331414
+      X=dMIN1(X,330.d0)
+      XK1=dSQRT(X)*dEXP(-X)*P
+      endif
+      RETURN
+      END
+
+
+
+
+c     ----------------------------------------------------------------------------------------------
+c     Added by RDW for inclusion of Baranov (2004) data
+      subroutine baranov(temp,wn,abcoef)
+      ! computes the dimer spectrum
+
+      implicit none
+      real*8 temp               !temperature
+      real*8 wn                 !wavenumber
+      real*8 abcoef             !absorption coefficient
+
+      integer nS,nT
+!      parameter(nS=1713)
+      parameter(nS=3000)
+      parameter(nT=9)
+
+      real*8 wn_arr(nS)
+      real*8 temp_arr(nT)
+      real*8 dim_arr(nS,nT)
+
+      integer pop
+      logical firstcall
+
+      save wn_arr, temp_arr, dim_arr
+
+      character*100 dt_file
+      integer strlen,ios
+      character*100 label
+      label='subroutine baranov'
+
+      dt_file='/san/home/rdword/kspectrum/line_data/CO2_dimer_data_NEW'
+!      dt_file='./data/CO2_dimer_data'
+      open(33,file=dt_file(1:strlen(dt_file)),
+     &     form='unformatted',
+     &     status='old',iostat=ios)
+      if (ios.ne.0) then        ! file not found
+         write(*,*) 'Error from ',label(1:strlen(label)),' :'
+         write(*,*) 'Data file could not be found:'
+         write(*,*) dt_file(1:strlen(dt_file))
+         stop
+      else
+         read(33) wn_arr
+         read(33) temp_arr
+         read(33) dim_arr
+      endif
+      close(33)
+
+
+!	print*,wn_arr(1)!
+!	print*,wn_arr(3000)
+!	print*,temp_arr
+!	print*,dim_arr(1500,1)
+!	stop
+
+
+      call bilinear(wn_arr,temp_arr,nS,nT,dim_arr,wn,temp,abcoef)
+
+ 111  continue
+      return
+      end
+c     ----------------------------------------------------------------------------------------------
+
+c     ----------------------------------------------------------------------------------------------
+c     Added by RDW for inclusion of Baranov (2004) data
+      subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x,y,f)
+
+      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_arr(nX)
+      real*8 y_arr(nY)
+      real*8 f2d_arr(nX,nY)
+      
+      integer strlen
+      character*100 label
+      label='subroutine bilinear'
+
+c     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
+c     print*,'x=',x,x_Arr(2),x_arr(nX-2)
+c     stop('Outside CIA dimer wavenumber range!\n')
+         goto 112
+      else
+         
+c     in the x (wavenumber) direction 1st
+         i=1
+ 10      if (i.lt.(nX+1)) then ! what passes for a 'while' loop in this sorry language...
+            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(*,*) 'Error from ',label(1:strlen(label)),' :'
+         write(*,*) 'Outside CIA dimer 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
+!         stop
+      else
+c     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)
+      
+c     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)
+      
+c     then in y-direction
+      f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
+      
+ 112  continue
+      return
+      end
+c     ----------------------------------------------------------------------------------------------
+      
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/compile_corrk
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/compile_corrk	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/compile_corrk	(revision 1512)
@@ -0,0 +1,11 @@
+#!/bin/bash
+#gfortran -c -g -fbounds-check generate_kmatrix.F90 strlen.for
+
+gfortran -c generate_kmatrix.F90 co2_cont.for strlen.for
+gfortran -o generate_kmatrix.exe generate_kmatrix.o co2_cont.o strlen.o
+
+gfortran -c make_composition_gen.for strlen.for
+gfortran -o make_composition_gen.exe make_composition_gen.o strlen.o
+
+rm -f *o
+rm -f *~
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/38x36/narrowbands_IR.in
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/38x36/narrowbands_IR.in	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/38x36/narrowbands_IR.in	(revision 1512)
@@ -0,0 +1,39 @@
+ 38
+  10.0    100.0
+  100.0   160.0
+  160.0   220.0
+  220.0   280.0
+  280.0   330.0
+  330.0   380.0
+  380.0   440.0
+  440.0   495.0
+  495.0   545.0
+  545.0   617.0
+  617.0   667.0
+  667.0   720.0
+  720.0   800.0
+  800.0   875.0
+  875.0   940.0
+  940.0   1000.0
+  1000.0  1065.0
+  1065.0  1108.0
+  1108.0  1200.0
+  1200.0  1275.0
+  1275.0  1350.0
+  1350.0  1450.0
+  1450.0  1550.0
+  1550.0  1650.0
+  1650.0  1750.0
+  1750.0  1850.0
+  1850.0  1950.0
+  1950.0  2150.0
+  2150.0  2250.0
+  2250.0  2450.0
+  2450.0  2650.0
+  2650.0  2850.0
+  2850.0  3000.0
+  3000.0  3200.0
+  3200.0  3400.0
+  3400.0  3600.0
+  3600.0  3800.0
+  3800.0  4000.0
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/38x36/narrowbands_VI.in
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/38x36/narrowbands_VI.in	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/38x36/narrowbands_VI.in	(revision 1512)
@@ -0,0 +1,37 @@
+           36
+    1950.000000000000         2050.000000000000     
+    2050.000000000000         2200.000000000000     
+    2200.000000000000         2397.000000000000     
+    2397.000000000000         2494.000000000000     
+    2494.000000000000         2796.000000000000     
+    2796.000000000000         3087.000000000000     
+    3087.000000000000         3425.000000000000     
+    3425.000000000000         3760.000000000000     
+    3760.000000000000         4030.000000000000     
+    4030.000000000000         4540.000000000000     
+    4540.000000000000         4950.000000000000     
+    4950.000000000000         5370.000000000000     
+    5370.000000000000         5925.000000000000     
+    5925.000000000000         6390.000000000000     
+    6390.000000000000         6990.000000000000     
+    6990.000000000000         7650.000000000000     
+    7650.000000000000         8315.000000000000     
+    8315.000000000000         8850.000000000000     
+    8850.000000000000         9350.000000000000     
+    9350.000000000000         9650.000000000000     
+    9650.000000000000         10400.00000000000     
+    10400.00000000000         11220.00000000000     
+    11220.00000000000         11870.00000000000     
+    11870.00000000000         12750.00000000000     
+    12750.00000000000         13300.00000000000     
+    13300.00000000000         14470.00000000000     
+    14470.00000000000         15000.00000000000     
+    15000.00000000000         15385.00000000000     
+    15385.00000000000         16667.00000000000     
+    16667.00000000000         18182.00000000000     
+    18182.00000000000         20000.00000000000     
+    20000.00000000000         22222.00000000000     
+    22222.00000000000         25000.00000000000     
+    25000.00000000000         28571.00000000000     
+    28571.00000000000         33333.00000000000     
+    33333.00000000000         35000.00000000000
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/Q.dat
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/Q.dat	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/Q.dat	(revision 1512)
@@ -0,0 +1,15 @@
+	4
+N2 
+CO2
+CH4
+H2O
+	7
+1.0e-7
+1.0e-6
+1.0e-5
+1.0e-4
+1.0e-3
+1.0e-2
+1.0e-1
+
+
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/T.dat
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/T.dat	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/T.dat	(revision 1512)
@@ -0,0 +1,7 @@
+	6
+100
+150
+200
+250
+300
+350
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/p.dat
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/p.dat	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/corrk_test/p.dat	(revision 1512)
@@ -0,0 +1,10 @@
+	9
+-3
+-2
+-1
+0
+1
+2
+3
+4
+5
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/formats.inc
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/formats.inc	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/formats.inc	(revision 1512)
@@ -0,0 +1,42 @@
+ 10   format(a)
+ 11   format(I1)
+ 12   format(I2)
+ 13   format(I3)
+ 14   format(I4)
+ 15   format(I5)
+ 16   format(I6)
+ 17   format(I7)
+ 18   format(I8)
+ 19   format(I9)
+ 110   format(I10)
+ 20   format(a23,I3,I3)
+ 21   format(a1)
+ 22   format(i3,a1,a10)
+ 30   format(a30,i3)
+ 31   format(a18,i3)
+ 32   format(a22,i3,a1,i3)
+ 33   format(a23,i3,a1,i3)
+ 34   format(a21,i3)
+ 35   format(a36,i3)
+ 36   format(a3,e17.9,a9,i4)
+   37   format(I1,I9,a100)
+   38   format(I9,a100)
+   39   format(a30,i3)
+   40  format(a,f9.3,a,f9.3,a)
+   41 format(a,f9.3,a)
+   42 format(a,I3,a,I3,a,I3,a)
+   43 format(a,I3,a)
+   44 format(I2,a)
+   45 format(I2,a,I2,a)
+   46 format(I2,a,I2,a,I2,a)
+   47 format(I2,a,I2,a,I2,a,I2,a)
+ 50   format(1000e17.9)
+ 51   format(I3,1000e17.9)
+ 52   format(1000e23.15)
+ 60   format(I2,I1,F12.6,2E10.3,2F5.4,F10.4,F4.2,F8.6,! &
+     &       4A15,6I1,6I2,A1,2F7.1)
+ 62   format(i12,F13.5,F14.4,i5,E14.6)
+ 63   format(I2)
+ 70   format(I2,I1,F12.6,1P2E10.3,0P2F5.4,F10.4,F4.2,F8.6,! &
+     &       2I3,2A9,3I1,3I2)
+ 80   format(I2,I1,F12.6)
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/generate_kmatrix.F90
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/generate_kmatrix.F90	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/generate_kmatrix.F90	(revision 1512)
@@ -0,0 +1,461 @@
+      program generate_kmatrix
+      implicit none
+
+!     -------------------------------------------------------------
+!     Purpose: to read hires spectra and produce kdistribution data
+!     Authour: Robin Wordsworth (2010)
+!     -------------------------------------------------------------
+
+      integer nmolec,mol
+      integer i,j,k,im
+      integer Nb
+
+      character*4 fname
+      character*32 fnum
+
+      integer iT, iP, iQ, Nlevs
+
+      double precision nuT, sigT, kT
+
+      integer Nq
+      parameter (Nq=16) ! we choose this here
+
+      integer Nmax
+      parameter (Nmax=1000) ! maximum number of levels
+
+      integer Nmol_max
+      parameter (Nmol_max=10) ! maximum number of radiatively active species
+
+      integer Nkmax
+      ! large array = slow calculation... this makes a huge difference
+      !parameter (Nkmax=50000) ! works for IR, 27 bands
+      !parameter (Nkmax=100000) ! works for IR, 27 bands
+      !parameter (Nkmax=140000) ! works for IR, 32 bands
+      !parameter (Nkmax=250000) ! works for VI, 36 bands
+      parameter (Nkmax=400000) ! works for IR, 8 bands
+      !parameter (Nkmax=400000) ! works for VI, 12 bands
+      !parameter (Nkmax=200000) ! works for IR, 30 bands
+
+      integer NbMAX
+      parameter (NbMAX=100)
+
+      integer Ngres
+      parameter (Ngres=200)
+
+      character*10 molec_names(1:Nmol_max)
+      double precision T(1:Nmax), P(1:Nmax), Q(1:Nmax)
+
+
+      double precision pi
+      parameter (pi=3.14159265)
+      double precision x
+      double precision y(1:Nq+2)
+
+      double precision nu_min(1:NbMAX) ! make allocatable?
+      double precision nu_max(1:NbMAX)
+
+      integer quad, band
+
+      double precision ktemp(1:Nq)
+      double precision w(1:Nq), gw(1:Nq)
+
+      double precision, allocatable :: kdist_temp(:)
+      double precision, allocatable :: kdist(:,:,:,:,:)
+
+      double precision nuband(1:Nkmax), kband(1:Nkmax), nuwband(1:Nkmax) ! better static I think
+      double precision nuwbtot
+
+      double precision gfine(1:Ngres)      
+      double precision klevs(1:Ngres)
+      double precision addline(1:Nkmax)
+
+      double precision kmax, kmin, lkmax, lkmin, dlogK, klev
+      integer ic, ib, ig, ik1, ik2
+
+      integer filerr
+      logical filexi
+
+      logical addCIA
+      integer iCO2
+      double precision k1,k2,kCIA,amg,Pp
+
+      real truemean, distmean, etot
+      
+
+      addCIA=.false.
+
+      print*,'Creating the GCM correlated-k data...'
+      print*,' '
+
+      !--------------------------------------------------------
+      ! load temperature, pressure, mixing ratio data
+      open(9,file='T.dat')
+      read(9,*) iT
+      do i=1,iT
+         read(9,*) T(i)
+      enddo
+      close(9)
+
+      open(10,file='p.dat')
+      read(10,*) iP
+      do i=1,iP
+         read(10,*) P(i)
+         P(i)=10.**P(i)  ! note we keep P in millibars for the conversion later
+      enddo
+      close(10)
+
+      open(11,file='Q.dat')
+      read(11,*) nmolec
+      do mol=1,nmolec
+         read(11,*) molec_names(mol)
+      enddo
+      read(11,*) iQ
+      do i=1,iQ
+         read(11,*) Q(i)
+      enddo
+      close(11)
+
+      Nlevs=iT*iP*iQ ! total number of layers
+
+      print*,'Temperature layers:  ',iT
+      print*,'Pressure layers:     ',iP
+      print*,'Mixing ratio layers: ',iQ
+      print*,'Total:               ',Nlevs
+
+
+      if(addCIA)then
+
+         if(nmolec.gt.2)then
+            print*,'I cant deal with this situation yet'
+            call abort
+         endif
+
+         iCO2=0
+         do mol=1,nmolec
+            if(molec_names(mol).eq.'CO2')then
+               iCO2=mol
+            endif
+         enddo
+
+         !if(iCO2.eq.1)then
+         print*,' '
+         print*,'--- Adding CO2 CIA ---'
+         !else
+         !   print*,'Warning, not CO2+H2O mixture, need CO2 mixing ratio'
+         !   call abort
+         !endif
+      endif
+
+      !--------------------------------------------------------
+      ! load spectral data
+      open(12,file='narrowbands.in')
+      read(12,*) Nb
+
+      do band=1,Nb
+         read(12,*) nu_min(band),nu_max(band)
+      enddo
+      close(12)
+
+      print*,' '
+      print*,'Number of spectral bands: ',Nb
+
+      print*,'Band limits:'
+      do band=1,Nb
+         print*,band,'-->',nu_min(band),nu_max(band),' cm^-1'
+      end do
+
+      !--------------------------------------------------------
+      ! create g-space data
+      do quad=0,Nq+1
+         x=dble(quad)/dble(Nq+2)
+         !y(quad+1)=sin(pi*x)*exp(-5*x)
+         !y(quad+1)=sin(pi*x)*exp(-9*x)  ! CO2_H2Ovar 38x36
+         y(quad+1)=sin(pi*x)*exp(-12*x)
+         !y(quad+1)=sin(pi*x)*exp(-15*x) ! CO2_CH4
+      enddo
+
+      do quad=1,Nq
+         w(quad)=y(quad+1)
+      enddo
+      w=w/sum(w)
+
+      print*,' '
+      print*,'Number of g-space intervals: ',Nq
+
+      print*,'g-space weighting:'
+      do quad=1,Nq
+         print*,quad,'-->',w(quad)
+      end do
+
+      ! write g-space data
+      open(14,file='g.dat')
+      write(14,*) Nq+1
+
+      do quad=1,Nq
+         write(14,*) w(quad)
+      enddo
+      write(14,*) 0.0
+      close(14)
+
+      ! cumulative sum for g-vector
+      gw(1)=w(1)
+      do quad=2,Nq
+         gw(quad)=gw(quad-1)+w(quad)
+      enddo
+
+      ! allocate 5D matrices
+      allocate(kdist_temp(Nq))
+      allocate(kdist(iT,iP,iQ,Nb,Nq+1))
+
+
+      nuT=0.0
+      do i=1,iQ
+         do j=1,iP
+
+
+            if(addCIA.and.(iCO2.eq.1))then 
+               Pp=P(j)*(1.0-Q(i))
+            else
+               Pp=P(j)*Q(i)
+            endif
+            !if(addCIA) Pp=P(j)*(1.0-Q(i))
+
+
+            do k=1,iT             
+               im = (i-1)*iP*iT + (j-1)*iT + k
+
+               !-----------------------------------------------
+               ! open kspectrum file
+               write(fnum,*) im
+               if(im<10)then
+                   fname='k00'//trim(adjustl(fnum))
+               elseif(im<100)then
+                   fname='k0'//trim(adjustl(fnum))
+               else
+                   fname='k'//trim(adjustl(fnum))
+               endif
+
+               ! check that the file exists
+               inquire(FILE='hires_spectrum/'//fname,EXIST=filexi)
+               if(.not.filexi) then
+                  write(*,*)'The file ',fname(1:LEN_TRIM(fname))
+                  write(*,*)'was not found in the folder hires_spectrum, exiting.'
+                  call abort
+               endif
+
+               print*,'Opening file ',fname(1:LEN_TRIM(fname))
+               open(17,file='hires_spectrum/'//fname)
+
+               nuT=0.0
+               do band=1,Nb
+               print*,'Band ---> ',band
+
+                  !--------------------------------------------
+                  ! write kdata to temporary band arrays
+                  ! and find max / min k-values in band
+                  do while (nuT.lt.nu_min(band))
+                      read(17,*,IOSTAT=filerr) nuT, sigT, kT
+                   ! ideally need to catch all errors here (ask Ehouarn)
+                      if(filerr.lt.0)then
+                         print*,    &
+                             'nu_min greater than highest value in file'
+                         print*,'IOSTAT=',filerr
+                         print*,'nu_min=',nu_min(band)
+                         print*,'nu_endoffile=',nuT
+                         call abort
+                      endif
+                  enddo
+
+                  ic=1
+                  kmin=kT
+                  kmax=kT
+                  nuband(:)=0
+                  kband(:)=0
+
+                  do while (nuT.lt.nu_max(band))
+
+                      nuband(ic)=nuT
+                      kband(ic)=kT
+
+                      ! insert CIA option here
+                      if(addCIA)then
+
+                         k1=0.0
+                         k2=0.0
+                         if(nuT.lt.500.0) call getspc(T(k),nuT,k1)
+                         if((nuT.gt.1000.0).and.(nuT.lt.1800.0)) call baranov(T(k),nuT,k2)
+                         ! cm^-1.amagat^2
+
+                         amg=273.15D+0/T(k)*(Pp/1013.25) ! amagats of CO2
+
+                         kCIA=(k1+k2)*1.0D+2*amg**2.0D+0 ! m^-1
+                         kband(ic)=kband(ic)+kCIA
+                   	 kT=kband(ic)
+
+                         !kband(ic)=kCIA
+			 !kT=kCIA
+			 !if(nuband(ic).gt.1250.0)then
+			 !	print*,'P=',P(i)
+			 !	print*,'T=',T(k)
+		         !	print*,'amg=',amg
+			 !	print*,'nu=',nuband(ic)
+                         !       print*,'k1=',k1
+                         !       print*,'k2=',k2
+                         !       print*,'kCIA=',kCIA
+		  	 !	call abort
+                         !endif
+
+                      endif
+
+                      if(kT.lt.kmin)then
+                        kmin=kT
+                      endif
+                      if(kT.gt.kmax)then ! bug corrected!
+                        kmax=kT
+                      endif
+
+                      read(17,*,IOSTAT=filerr) nuT, sigT, kT
+                      if(filerr.lt.0)then
+                         print*,    &
+                             'nu_max greater than highest value in file'
+                         print*,'IOSTAT=',filerr
+                         print*,'nu_max=',nu_max(band)
+                         print*,'nu_endoffile=',nuT
+                         call abort
+                      endif
+
+
+                      ic=ic+1
+                      
+                      if(ic.gt.Nkmax)then
+                         print*,'Spectral resolution too low for these bands, exiting!'
+                         print*,'Increase Nkmax in generate_kmatrix.F90.'
+                         call abort
+                      endif
+
+                  enddo
+                  nuband(ic)=nuT
+                  kband(ic)=kT
+
+                  if(kmax.le.0.0)then
+                  !if(kmax.le.1.0e-15)then
+                     print*,'kmax=',kmax
+                     print*,'Setting values to zero in this band.'
+                     kdist(k,j,i,band,1:(Nq+1))=0.0
+                  else
+
+                  if(kmin.lt.1e-30)then
+                     kmin=1.0e-30
+                  endif
+
+                  !--------------------------------------------
+                  ! calculate delta nu values
+                  nuwband(:)=0.0
+                  do ib=2,(ic-1)
+                      nuwband(ib)=(nuband(ib+1)-nuband(ib-1))/2
+                  enddo
+                  nuwband(1)=nuwband(2)
+                  nuwband(ic)=nuwband(ic-1) ! causes tiny boundary error; whatever
+                  nuwbtot=sum(nuwband)
+
+                  lkmax=log10(kmax)
+                  lkmin=log10(kmin)
+                  
+                  dlogK=(lkmax-lkmin)/(Ngres-1)                 
+
+                  !--------------------------------------------
+                  ! calculate the g function
+                  gfine(:)=0.0
+                  do ig=1,Ngres
+                    klev=lkmin + dble(ig)*dlogK
+                    klev=10**klev
+                    klevs(ig)=klev
+
+                    addline(:)=0
+                    do ib=1,ic
+                        if(kband(ib).lt.klev)then
+                           addline(ib)=1.0
+                        endif
+                    enddo
+                    gfine(ig)=sum(addline(:)*nuwband(:))/nuwbtot
+
+                  enddo
+                  print*,'gnorm    = ',gfine(Ngres)
+
+                  !--------------------------------------------
+                  ! invert the g function
+                  ik1=1
+                  ik2=1
+                  kdist_temp(:)=0.0
+                  do quad=1,Nq
+                    
+                    ik2=ik1
+
+                    gfind: do while (gfine(ik2).lt.gw(quad))
+                       ik2=ik2+1
+                       if(ik2.gt.Ngres)then
+                          print*,'g-inversion overflow...'
+                          ik2=Ngres
+                          exit gfind
+                       endif
+                    enddo gfind
+
+                    !if(ik2.gt.Ngres)then ! avoid out-of-bounds on final step
+                    !   ik2=Ngres
+                    !endif
+
+
+                     if(.not.(ik1.ge.ik2))then ! this could probably be improved
+                        !kdist_temp(quad)=sum(klevs(ik1:ik2))/(ik2-ik1) ! a bug was here!
+                        kdist_temp(quad)=sum(klevs(ik1:ik2))/(ik2+1-ik1)
+
+
+                     elseif(quad.ne.1)then
+                        kdist_temp(quad)=kdist_temp(quad-1)
+                     else
+                        kdist_temp(quad)=kmin
+                     endif
+
+                     !-----------------------------------------
+                     ! assign to 5D matrix and convert to cm^2 / molecule
+                     kdist(k,j,i,band,quad)=1.3807e-21 *               &
+                                kdist_temp(quad) * T(k)/P(j)
+
+                     ik1=ik2
+                  enddo ! g-space
+                  kdist(k,j,i,band,Nq+1)=0.0
+
+                  truemean = sum(kband(1:ic)*nuwband(1:ic))/sum(nuwband(1:ic))
+                  distmean = sum(kdist_temp(:)*w(:))/sum(w)
+                  etot     = 100*abs(truemean-distmean)/(distmean+truemean)
+                  print*,'truemean = ',truemean
+                  print*,'distmean = ',distmean
+                  !print*,'kdisttemp = ',kdist_temp(Nq)
+                  !print*,'kdistmax  = ',kmax
+                  print*,'etot (%) = ',etot
+                  !print*,'w  = ',w
+
+                  endif ! if kmax .le. 0.0 
+
+               enddo ! bands
+
+            enddo
+         enddo
+      enddo
+
+      !--------------------------------------------------------
+      ! save data and clean up
+      close(17)
+      open(14,file='corrk_gcm.dat',form='formatted')
+      write(14,*) kdist
+      close(14)
+      print*,' '
+      print*,'Correlated-k data saved in corrk_gcm.dat'
+      print*,' '
+      print*,'Now you must change the values of L_NTREF etc. in radinc_h.F90 (if necessary)'
+      print*,'and the variable "corrkdata" in callphys.def'
+
+      deallocate(kdist_temp)
+      deallocate(kdist)
+
+    end program generate_kmatrix
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/make_composition_gen.for
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/make_composition_gen.for	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/make_composition_gen.for	(revision 1512)
@@ -0,0 +1,161 @@
+      program make_composition_gen
+      implicit none
+
+!     -------------------------------------------------------------   
+!     Purpose: to create "composition.in" file readable by kspectrum
+!     Authour: Robin Wordsworth (2010)
+!     -------------------------------------------------------------
+
+      include 'max.inc'
+      include 'formats.inc'
+      integer strlen,nmolec,nlev_mf,mol,imf,ind,indf
+      integer i,j,k,im,m
+      integer nb,na,ne
+      double precision alt(1:Nmax)
+      double precision pres(1:Nmax)
+      double precision temp(1:Nmax)
+
+      double precision P_mf(1:Nmax)
+      double precision x(1:Nmax,1:Nmol_max)
+      double precision xlev(1:Nmol_max)
+      double precision xconst(1:Nmol_max)
+      character*100 composition_file,model_file
+      character*100 planet_descriptor
+      character*200 string
+      character*20 s
+      character*10 molec_names(1:Nmol_max),molname
+      character*1 spch
+
+      integer iP, iT, iQ
+      double precision P(1:Nmax), T(1:Nmax), Q(1:Nmax)
+
+      spch=' '
+
+      print*,'Name of atmosphere / planet:'
+      read*,planet_descriptor
+
+      composition_file='./composition.in'
+
+!     ! nmolec=2 ! total number of molecules
+
+      ! load temperature, pressure, mixing ratio data
+      open(9,file='T.dat')
+      read(9,*) iT
+      do i=1,iT
+         read(9,*) T(i)
+      enddo
+      close(9)
+
+      open(10,file='p.dat')
+      read(10,*) iP
+      do i=1,iP
+         read(10,*) P(i)
+         P(i)=10.**P(i)/1013.25
+      enddo
+      close(10)
+
+      open(11,file='Q.dat')
+      read(11,*) nmolec
+
+      print*,'nmolec=',nmolec
+
+      do mol=1,nmolec
+         read(11,*) molec_names(mol)
+      enddo
+      read(11,*) iQ
+      do i=1,iQ
+         read(11,*) Q(i)
+      enddo
+      close(11)
+
+      m=iP*iT*iQ ! total number of layers
+
+      print*,'Temperature layers:  ',iT
+      print*,'Pressure layers:     ',iP
+      print*,'Mixing ratio layers: ',iQ
+      print*,'Total:               ',m
+
+      do mol=1,nmolec-1
+         print*,'Please enter vmr of ',molec_names(mol)
+         read(*,*) xconst(mol)
+      end do
+
+     
+      do i=1,iQ
+         do j=1,iP              
+            do k=1,iT
+
+               im = (i-1)*iP*iT + (j-1)*iT + k
+
+               alt(im)=0.D0
+               pres(im)=P(j)
+               temp(im)=T(k)
+
+               do mol=1,nmolec-1
+                  x(im,mol)=xconst(mol)*(1.0-Q(i))
+               enddo
+               x(im,nmolec)=Q(i)
+
+            enddo
+         enddo
+      enddo
+
+      open(12,file=composition_file(1:strlen(composition_file)))
+      string='Atmospheric composition input '
+      string=string(1:strlen(string))//' data file for planet:'
+      string=string(1:strlen(string))//' '//
+
+     &     planet_descriptor(1:strlen(planet_descriptor))
+      write(12,'(a)') string(1:strlen(string))
+      write(12,39) 'Number of atmospheric levels: ',m
+      write(12,34) 'Number of molecules: ',nmolec
+      write(12,*) 
+      string='      z (km)     /    P (atm)     /     T (K) *'
+
+
+      do mol=1,nmolec
+
+         molname=molec_names(mol)
+
+         nb=5
+         if (strlen(molname).eq.1) then
+            na=7
+         else if (strlen(molname).eq.2) then
+            na=6
+         else if (strlen(molname).eq.3) then
+            na=5
+         endif
+         s='*'
+         do ne=1,nb
+            s=s(1:strlen(s)-1)//spch//'*'
+         enddo
+         s=s(1:strlen(s)-1)//'/*'
+         do ne=1,na
+            s=s(1:strlen(s)-1)//spch//'*'
+         enddo
+         
+         string=string(1:strlen(string)-1)
+     &        //s(1:strlen(s)-1)
+     &        //'x['
+     &        //molname(1:strlen(molname))
+     &        //']*'
+      enddo
+      write(12,*) string(1:strlen(string)-1)
+
+      do im=1,m
+
+         do mol=1,nmolec
+            xlev(mol)=x(im,mol)
+         enddo                  ! mol
+
+         write(12,50) alt(im),pres(im),temp(im),(xlev(mol),mol=1,nmolec)
+
+      enddo
+
+      write(12,*) 
+      close(12)
+
+      write(*,*) 'Output file successfully generated:'
+      write(*,*) composition_file(1:strlen(composition_file))
+      
+      end
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/max.inc
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/max.inc	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/max.inc	(revision 1512)
@@ -0,0 +1,23 @@
+          integer Nline_mx,Nmol_max,Niso_max,Nbmx,Nkmx,Nmax,Nlpnb_mx
+	    integer Nxtab_mx,Nacc_mx,Nmax_tab,Nytab_mx,Nktab_mx
+	    integer nfiles_mx,Nproc_mx,Nchunk_mx,Nc_mx,Nmax_vt,Nppc_max
+	    integer Nti_mx
+          parameter(Nline_mx=2000000)
+          parameter(Nmol_max=42)
+          parameter(Niso_max=15)
+	  parameter(Nbmx=10000)
+	  parameter(Nkmx=100000)
+	  parameter(Nmax=2000)
+	  parameter(Nlpnb_mx=500000)
+	  parameter(Nxtab_mx=100000)
+	  parameter(Nacc_mx=1000)
+	  parameter(Nmax_tab=1000000)
+	    parameter(Nytab_mx=210)
+	    parameter(Nktab_mx=1000000)
+	    parameter(nfiles_mx=100)
+	    parameter(Nproc_mx=1000)
+	    parameter(Nchunk_mx=1000)
+	    parameter(Nc_mx=9)
+	    parameter(Nmax_vt=1000000)
+	    parameter(Nppc_max=10000)
+	    parameter(Nti_mx=100000)
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/postkspectrum
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/postkspectrum	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/postkspectrum	(revision 1512)
@@ -0,0 +1,48 @@
+#!/bin/bash
+
+DWORK_DIR=$RAD/CORRELATEDK/avec_robin_data
+KSPEC_DIR=$KSPECTRUM/kspec_run
+#BANDS_DIR=32x36
+BANDS_DIR=38x36
+
+### if hires_spectrum don't exist then create it, copy data from kspectrum
+DIRCHECK=hires_spectrum
+cd $DWORK_DIR/
+if [ -d $DIRCHECK ]; then
+    echo -en "Directory hires_spectrum already exists.\n"
+else
+    echo -en "Copying kspectrum data to its new home...\n"
+    cd $DWORK_DIR
+    mkdir hires_spectrum
+
+    mv $KSPEC_DIR/results/k* hires_spectrum/
+    cp $KSPEC_DIR/results/calculation_info.txt hires_spectrum/
+    cp $KSPEC_DIR/*in hires_spectrum/
+    cp $KSPEC_DIR/data/composition.in hires_spectrum/
+
+    # removing contents of 'optimizations' folder
+    rm -rf $KSPEC_DIR/optimizations
+    mkdir $KSPEC_DIR/optimizations
+
+    chmod 444 hires_spectrum/*
+fi
+
+cp ~/utilities/run_kmatrix .
+
+echo -en "\nNow you must submit the parallel job (>> llsubmit run_kmatrix)\n"
+
+exit 0
+
+#cp ~/utilities/generate_kmatrix.exe .
+### make kmatrix in the infrared and save it in a GCM-friendly format
+#cp $BANDS_DIR/narrowbands_IR.in narrowbands.in
+#generate_kmatrix.exe
+#mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_IR.dat 
+#rm narrowbands.in
+### make kmatrix in the visible and save it in a GCM-friendly format
+#cp $BANDS_DIR/narrowbands_VI.in narrowbands.in
+#generate_kmatrix.exe
+#mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_VI.dat 
+#rm narrowbands.in
+#rm generate_kmatrix.exe
+#exit 0
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/prekspectrum
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/prekspectrum	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/prekspectrum	(revision 1512)
@@ -0,0 +1,28 @@
+#!/bin/bash
+
+#DWORK_DIR=$RAD/CORRELATEDK/premier_essai
+
+DWORK_DIR=$PWD
+KSPEC_DIR=$KSPECTRUM/kspec_run
+
+cd $DWORK_DIR
+
+# test for presence of k* in kspec results directory
+FILECHECK=$KSPEC_DIR/results/k001
+if [ -f $FILECHECK ]; then
+    echo -en "There is already k-data in the kspec directory "
+    echo -en $KSPEC_DIR
+    echo -en ", delete or move it before starting!\n"
+    exit 1
+fi
+
+make_composition_gen.exe
+cp composition.in $KSPEC_DIR/data/
+
+echo -en "kspectrum working directory is:\n"
+echo -en $KSPEC_DIR
+echo -en "\nNow you must submit the parallel job (>> llsubmit run_kspec)\n"
+
+cd $KSPEC_DIR
+
+exit 0
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/run_kmatrix
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/run_kmatrix	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/run_kmatrix	(revision 1512)
@@ -0,0 +1,28 @@
+# @ job_name=kmatrix
+# @ output = $(job_name).$(jobid)
+# @ error = $(output)
+# @ job_type = serial
+# @ class = B
+# @ queue
+
+BANDS_DIR=38x36
+
+cp $MOD/LMDZ.GENERIC/utilities/generate_kmatrix.exe .
+
+## make kmatrix in the infrared and save it in a GCM-friendly format
+cp $BANDS_DIR/narrowbands_IR.in narrowbands.in
+generate_kmatrix.exe > log.out
+mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_IR.dat 
+mv log.out $BANDS_DIR/logIR.out 
+rm narrowbands.in
+
+## make kmatrix in the visible and save it in a GCM-friendly format
+cp $BANDS_DIR/narrowbands_VI.in narrowbands.in
+generate_kmatrix.exe > log.out
+mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_VI.dat 
+mv log.out $BANDS_DIR/logVI.out 
+rm narrowbands.in
+
+rm generate_kmatrix.exe
+
+exit 0
Index: trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/strlen.for
===================================================================
--- trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/strlen.for	(revision 1512)
+++ trunk/LMDZ.GENERIC/utilities/correlatedk_utilities/strlen.for	(revision 1512)
@@ -0,0 +1,10 @@
+      integer function strlen(st)
+      integer i
+      character st*(*)
+      i=len(st)
+      do while (st(i:i).eq.' ')
+         i=i-1
+      enddo
+      strlen=i
+      return
+      end
Index: trunk/LMDZ.GENERIC/utilities/formats.inc
===================================================================
--- trunk/LMDZ.GENERIC/utilities/formats.inc	(revision 1507)
+++ 	(revision )
@@ -1,42 +1,0 @@
- 10   format(a)
- 11   format(I1)
- 12   format(I2)
- 13   format(I3)
- 14   format(I4)
- 15   format(I5)
- 16   format(I6)
- 17   format(I7)
- 18   format(I8)
- 19   format(I9)
- 110   format(I10)
- 20   format(a23,I3,I3)
- 21   format(a1)
- 22   format(i3,a1,a10)
- 30   format(a30,i3)
- 31   format(a18,i3)
- 32   format(a22,i3,a1,i3)
- 33   format(a23,i3,a1,i3)
- 34   format(a21,i3)
- 35   format(a36,i3)
- 36   format(a3,e17.9,a9,i4)
-   37   format(I1,I9,a100)
-   38   format(I9,a100)
-   39   format(a30,i3)
-   40  format(a,f9.3,a,f9.3,a)
-   41 format(a,f9.3,a)
-   42 format(a,I3,a,I3,a,I3,a)
-   43 format(a,I3,a)
-   44 format(I2,a)
-   45 format(I2,a,I2,a)
-   46 format(I2,a,I2,a,I2,a)
-   47 format(I2,a,I2,a,I2,a,I2,a)
- 50   format(1000e17.9)
- 51   format(I3,1000e17.9)
- 52   format(1000e23.15)
- 60   format(I2,I1,F12.6,2E10.3,2F5.4,F10.4,F4.2,F8.6,! &
-     &       4A15,6I1,6I2,A1,2F7.1)
- 62   format(i12,F13.5,F14.4,i5,E14.6)
- 63   format(I2)
- 70   format(I2,I1,F12.6,1P2E10.3,0P2F5.4,F10.4,F4.2,F8.6,! &
-     &       2I3,2A9,3I1,3I2)
- 80   format(I2,I1,F12.6)
Index: trunk/LMDZ.GENERIC/utilities/generate_kmatrix.F90
===================================================================
--- trunk/LMDZ.GENERIC/utilities/generate_kmatrix.F90	(revision 1507)
+++ 	(revision )
@@ -1,461 +1,0 @@
-      program generate_kmatrix
-      implicit none
-
-!     -------------------------------------------------------------
-!     Purpose: to read hires spectra and produce kdistribution data
-!     Authour: Robin Wordsworth (2010)
-!     -------------------------------------------------------------
-
-      integer nmolec,mol
-      integer i,j,k,im
-      integer Nb
-
-      character*4 fname
-      character*32 fnum
-
-      integer iT, iP, iQ, Nlevs
-
-      double precision nuT, sigT, kT
-
-      integer Nq
-      parameter (Nq=16) ! we choose this here
-
-      integer Nmax
-      parameter (Nmax=1000) ! maximum number of levels
-
-      integer Nmol_max
-      parameter (Nmol_max=10) ! maximum number of radiatively active species
-
-      integer Nkmax
-      ! large array = slow calculation... this makes a huge difference
-      !parameter (Nkmax=50000) ! works for IR, 27 bands
-      !parameter (Nkmax=100000) ! works for IR, 27 bands
-      !parameter (Nkmax=140000) ! works for IR, 32 bands
-      !parameter (Nkmax=250000) ! works for VI, 36 bands
-      parameter (Nkmax=400000) ! works for IR, 8 bands
-      !parameter (Nkmax=400000) ! works for VI, 12 bands
-      !parameter (Nkmax=200000) ! works for IR, 30 bands
-
-      integer NbMAX
-      parameter (NbMAX=100)
-
-      integer Ngres
-      parameter (Ngres=200)
-
-      character*10 molec_names(1:Nmol_max)
-      double precision T(1:Nmax), P(1:Nmax), Q(1:Nmax)
-
-
-      double precision pi
-      parameter (pi=3.14159265)
-      double precision x
-      double precision y(1:Nq+2)
-
-      double precision nu_min(1:NbMAX) ! make allocatable?
-      double precision nu_max(1:NbMAX)
-
-      integer quad, band
-
-      double precision ktemp(1:Nq)
-      double precision w(1:Nq), gw(1:Nq)
-
-      double precision, allocatable :: kdist_temp(:)
-      double precision, allocatable :: kdist(:,:,:,:,:)
-
-      double precision nuband(1:Nkmax), kband(1:Nkmax), nuwband(1:Nkmax) ! better static I think
-      double precision nuwbtot
-
-      double precision gfine(1:Ngres)      
-      double precision klevs(1:Ngres)
-      double precision addline(1:Nkmax)
-
-      double precision kmax, kmin, lkmax, lkmin, dlogK, klev
-      integer ic, ib, ig, ik1, ik2
-
-      integer filerr
-      logical filexi
-
-      logical addCIA
-      integer iCO2
-      double precision k1,k2,kCIA,amg,Pp
-
-      real truemean, distmean, etot
-      
-
-      addCIA=.false.
-
-      print*,'Creating the GCM correlated-k data...'
-      print*,' '
-
-      !--------------------------------------------------------
-      ! load temperature, pressure, mixing ratio data
-      open(9,file='T.dat')
-      read(9,*) iT
-      do i=1,iT
-         read(9,*) T(i)
-      enddo
-      close(9)
-
-      open(10,file='p.dat')
-      read(10,*) iP
-      do i=1,iP
-         read(10,*) P(i)
-         P(i)=10.**P(i)  ! note we keep P in millibars for the conversion later
-      enddo
-      close(10)
-
-      open(11,file='Q.dat')
-      read(11,*) nmolec
-      do mol=1,nmolec
-         read(11,*) molec_names(mol)
-      enddo
-      read(11,*) iQ
-      do i=1,iQ
-         read(11,*) Q(i)
-      enddo
-      close(11)
-
-      Nlevs=iT*iP*iQ ! total number of layers
-
-      print*,'Temperature layers:  ',iT
-      print*,'Pressure layers:     ',iP
-      print*,'Mixing ratio layers: ',iQ
-      print*,'Total:               ',Nlevs
-
-
-      if(addCIA)then
-
-         if(nmolec.gt.2)then
-            print*,'I cant deal with this situation yet'
-            call abort
-         endif
-
-         iCO2=0
-         do mol=1,nmolec
-            if(molec_names(mol).eq.'CO2')then
-               iCO2=mol
-            endif
-         enddo
-
-         !if(iCO2.eq.1)then
-         print*,' '
-         print*,'--- Adding CO2 CIA ---'
-         !else
-         !   print*,'Warning, not CO2+H2O mixture, need CO2 mixing ratio'
-         !   call abort
-         !endif
-      endif
-
-      !--------------------------------------------------------
-      ! load spectral data
-      open(12,file='narrowbands.in')
-      read(12,*) Nb
-
-      do band=1,Nb
-         read(12,*) nu_min(band),nu_max(band)
-      enddo
-      close(12)
-
-      print*,' '
-      print*,'Number of spectral bands: ',Nb
-
-      print*,'Band limits:'
-      do band=1,Nb
-         print*,band,'-->',nu_min(band),nu_max(band),' cm^-1'
-      end do
-
-      !--------------------------------------------------------
-      ! create g-space data
-      do quad=0,Nq+1
-         x=dble(quad)/dble(Nq+2)
-         !y(quad+1)=sin(pi*x)*exp(-5*x)
-         !y(quad+1)=sin(pi*x)*exp(-9*x)  ! CO2_H2Ovar 38x36
-         y(quad+1)=sin(pi*x)*exp(-12*x)
-         !y(quad+1)=sin(pi*x)*exp(-15*x) ! CO2_CH4
-      enddo
-
-      do quad=1,Nq
-         w(quad)=y(quad+1)
-      enddo
-      w=w/sum(w)
-
-      print*,' '
-      print*,'Number of g-space intervals: ',Nq
-
-      print*,'g-space weighting:'
-      do quad=1,Nq
-         print*,quad,'-->',w(quad)
-      end do
-
-      ! write g-space data
-      open(14,file='g.dat')
-      write(14,*) Nq+1
-
-      do quad=1,Nq
-         write(14,*) w(quad)
-      enddo
-      write(14,*) 0.0
-      close(14)
-
-      ! cumulative sum for g-vector
-      gw(1)=w(1)
-      do quad=2,Nq
-         gw(quad)=gw(quad-1)+w(quad)
-      enddo
-
-      ! allocate 5D matrices
-      allocate(kdist_temp(Nq))
-      allocate(kdist(iT,iP,iQ,Nb,Nq+1))
-
-
-      nuT=0.0
-      do i=1,iQ
-         do j=1,iP
-
-
-            if(addCIA.and.(iCO2.eq.1))then 
-               Pp=P(j)*(1.0-Q(i))
-            else
-               Pp=P(j)*Q(i)
-            endif
-            !if(addCIA) Pp=P(j)*(1.0-Q(i))
-
-
-            do k=1,iT             
-               im = (i-1)*iP*iT + (j-1)*iT + k
-
-               !-----------------------------------------------
-               ! open kspectrum file
-               write(fnum,*) im
-               if(im<10)then
-                   fname='k00'//trim(adjustl(fnum))
-               elseif(im<100)then
-                   fname='k0'//trim(adjustl(fnum))
-               else
-                   fname='k'//trim(adjustl(fnum))
-               endif
-
-               ! check that the file exists
-               inquire(FILE='hires_spectrum/'//fname,EXIST=filexi)
-               if(.not.filexi) then
-                  write(*,*)'The file ',fname(1:LEN_TRIM(fname))
-                  write(*,*)'was not found in the folder hires_spectrum, exiting.'
-                  call abort
-               endif
-
-               print*,'Opening file ',fname(1:LEN_TRIM(fname))
-               open(17,file='hires_spectrum/'//fname)
-
-               nuT=0.0
-               do band=1,Nb
-               print*,'Band ---> ',band
-
-                  !--------------------------------------------
-                  ! write kdata to temporary band arrays
-                  ! and find max / min k-values in band
-                  do while (nuT.lt.nu_min(band))
-                      read(17,*,IOSTAT=filerr) nuT, sigT, kT
-                   ! ideally need to catch all errors here (ask Ehouarn)
-                      if(filerr.lt.0)then
-                         print*,    &
-                             'nu_min greater than highest value in file'
-                         print*,'IOSTAT=',filerr
-                         print*,'nu_min=',nu_min(band)
-                         print*,'nu_endoffile=',nuT
-                         call abort
-                      endif
-                  enddo
-
-                  ic=1
-                  kmin=kT
-                  kmax=kT
-                  nuband(:)=0
-                  kband(:)=0
-
-                  do while (nuT.lt.nu_max(band))
-
-                      nuband(ic)=nuT
-                      kband(ic)=kT
-
-                      ! insert CIA option here
-                      if(addCIA)then
-
-                         k1=0.0
-                         k2=0.0
-                         if(nuT.lt.500.0) call getspc(T(k),nuT,k1)
-                         if((nuT.gt.1000.0).and.(nuT.lt.1800.0)) call baranov(T(k),nuT,k2)
-                         ! cm^-1.amagat^2
-
-                         amg=273.15D+0/T(k)*(Pp/1013.25) ! amagats of CO2
-
-                         kCIA=(k1+k2)*1.0D+2*amg**2.0D+0 ! m^-1
-                         kband(ic)=kband(ic)+kCIA
-                   	 kT=kband(ic)
-
-                         !kband(ic)=kCIA
-			 !kT=kCIA
-			 !if(nuband(ic).gt.1250.0)then
-			 !	print*,'P=',P(i)
-			 !	print*,'T=',T(k)
-		         !	print*,'amg=',amg
-			 !	print*,'nu=',nuband(ic)
-                         !       print*,'k1=',k1
-                         !       print*,'k2=',k2
-                         !       print*,'kCIA=',kCIA
-		  	 !	call abort
-                         !endif
-
-                      endif
-
-                      if(kT.lt.kmin)then
-                        kmin=kT
-                      endif
-                      if(kT.gt.kmax)then ! bug corrected!
-                        kmax=kT
-                      endif
-
-                      read(17,*,IOSTAT=filerr) nuT, sigT, kT
-                      if(filerr.lt.0)then
-                         print*,    &
-                             'nu_max greater than highest value in file'
-                         print*,'IOSTAT=',filerr
-                         print*,'nu_max=',nu_max(band)
-                         print*,'nu_endoffile=',nuT
-                         call abort
-                      endif
-
-
-                      ic=ic+1
-                      
-                      if(ic.gt.Nkmax)then
-                         print*,'Spectral resolution too low for these bands, exiting!'
-                         print*,'Increase Nkmax in generate_kmatrix.F90.'
-                         call abort
-                      endif
-
-                  enddo
-                  nuband(ic)=nuT
-                  kband(ic)=kT
-
-                  if(kmax.le.0.0)then
-                  !if(kmax.le.1.0e-15)then
-                     print*,'kmax=',kmax
-                     print*,'Setting values to zero in this band.'
-                     kdist(k,j,i,band,1:(Nq+1))=0.0
-                  else
-
-                  if(kmin.lt.1e-30)then
-                     kmin=1.0e-30
-                  endif
-
-                  !--------------------------------------------
-                  ! calculate delta nu values
-                  nuwband(:)=0.0
-                  do ib=2,(ic-1)
-                      nuwband(ib)=(nuband(ib+1)-nuband(ib-1))/2
-                  enddo
-                  nuwband(1)=nuwband(2)
-                  nuwband(ic)=nuwband(ic-1) ! causes tiny boundary error; whatever
-                  nuwbtot=sum(nuwband)
-
-                  lkmax=log10(kmax)
-                  lkmin=log10(kmin)
-                  
-                  dlogK=(lkmax-lkmin)/(Ngres-1)                 
-
-                  !--------------------------------------------
-                  ! calculate the g function
-                  gfine(:)=0.0
-                  do ig=1,Ngres
-                    klev=lkmin + dble(ig)*dlogK
-                    klev=10**klev
-                    klevs(ig)=klev
-
-                    addline(:)=0
-                    do ib=1,ic
-                        if(kband(ib).lt.klev)then
-                           addline(ib)=1.0
-                        endif
-                    enddo
-                    gfine(ig)=sum(addline(:)*nuwband(:))/nuwbtot
-
-                  enddo
-                  print*,'gnorm    = ',gfine(Ngres)
-
-                  !--------------------------------------------
-                  ! invert the g function
-                  ik1=1
-                  ik2=1
-                  kdist_temp(:)=0.0
-                  do quad=1,Nq
-                    
-                    ik2=ik1
-
-                    gfind: do while (gfine(ik2).lt.gw(quad))
-                       ik2=ik2+1
-                       if(ik2.gt.Ngres)then
-                          print*,'g-inversion overflow...'
-                          ik2=Ngres
-                          exit gfind
-                       endif
-                    enddo gfind
-
-                    !if(ik2.gt.Ngres)then ! avoid out-of-bounds on final step
-                    !   ik2=Ngres
-                    !endif
-
-
-                     if(.not.(ik1.ge.ik2))then ! this could probably be improved
-                        !kdist_temp(quad)=sum(klevs(ik1:ik2))/(ik2-ik1) ! a bug was here!
-                        kdist_temp(quad)=sum(klevs(ik1:ik2))/(ik2+1-ik1)
-
-
-                     elseif(quad.ne.1)then
-                        kdist_temp(quad)=kdist_temp(quad-1)
-                     else
-                        kdist_temp(quad)=kmin
-                     endif
-
-                     !-----------------------------------------
-                     ! assign to 5D matrix and convert to cm^2 / molecule
-                     kdist(k,j,i,band,quad)=1.3807e-21 *               &
-                                kdist_temp(quad) * T(k)/P(j)
-
-                     ik1=ik2
-                  enddo ! g-space
-                  kdist(k,j,i,band,Nq+1)=0.0
-
-                  truemean = sum(kband(1:ic)*nuwband(1:ic))/sum(nuwband(1:ic))
-                  distmean = sum(kdist_temp(:)*w(:))/sum(w)
-                  etot     = 100*abs(truemean-distmean)/(distmean+truemean)
-                  print*,'truemean = ',truemean
-                  print*,'distmean = ',distmean
-                  !print*,'kdisttemp = ',kdist_temp(Nq)
-                  !print*,'kdistmax  = ',kmax
-                  print*,'etot (%) = ',etot
-                  !print*,'w  = ',w
-
-                  endif ! if kmax .le. 0.0 
-
-               enddo ! bands
-
-            enddo
-         enddo
-      enddo
-
-      !--------------------------------------------------------
-      ! save data and clean up
-      close(17)
-      open(14,file='corrk_gcm.dat',form='formatted')
-      write(14,*) kdist
-      close(14)
-      print*,' '
-      print*,'Correlated-k data saved in corrk_gcm.dat'
-      print*,' '
-      print*,'Now you must change the values of L_NTREF etc. in radinc_h.F90 (if necessary)'
-      print*,'and the variable "corrkdata" in callphys.def'
-
-      deallocate(kdist_temp)
-      deallocate(kdist)
-
-    end program generate_kmatrix
Index: trunk/LMDZ.GENERIC/utilities/make_composition_gen.for
===================================================================
--- trunk/LMDZ.GENERIC/utilities/make_composition_gen.for	(revision 1507)
+++ 	(revision )
@@ -1,161 +1,0 @@
-      program make_composition_gen
-      implicit none
-
-!     -------------------------------------------------------------   
-!     Purpose: to create "composition.in" file readable by kspectrum
-!     Authour: Robin Wordsworth (2010)
-!     -------------------------------------------------------------
-
-      include 'max.inc'
-      include 'formats.inc'
-      integer strlen,nmolec,nlev_mf,mol,imf,ind,indf
-      integer i,j,k,im,m
-      integer nb,na,ne
-      double precision alt(1:Nmax)
-      double precision pres(1:Nmax)
-      double precision temp(1:Nmax)
-
-      double precision P_mf(1:Nmax)
-      double precision x(1:Nmax,1:Nmol_max)
-      double precision xlev(1:Nmol_max)
-      double precision xconst(1:Nmol_max)
-      character*100 composition_file,model_file
-      character*100 planet_descriptor
-      character*200 string
-      character*20 s
-      character*10 molec_names(1:Nmol_max),molname
-      character*1 spch
-
-      integer iP, iT, iQ
-      double precision P(1:Nmax), T(1:Nmax), Q(1:Nmax)
-
-      spch=' '
-
-      print*,'Name of atmosphere / planet:'
-      read*,planet_descriptor
-
-      composition_file='./composition.in'
-
-!     ! nmolec=2 ! total number of molecules
-
-      ! load temperature, pressure, mixing ratio data
-      open(9,file='T.dat')
-      read(9,*) iT
-      do i=1,iT
-         read(9,*) T(i)
-      enddo
-      close(9)
-
-      open(10,file='p.dat')
-      read(10,*) iP
-      do i=1,iP
-         read(10,*) P(i)
-         P(i)=10.**P(i)/1013.25
-      enddo
-      close(10)
-
-      open(11,file='Q.dat')
-      read(11,*) nmolec
-
-      print*,'nmolec=',nmolec
-
-      do mol=1,nmolec
-         read(11,*) molec_names(mol)
-      enddo
-      read(11,*) iQ
-      do i=1,iQ
-         read(11,*) Q(i)
-      enddo
-      close(11)
-
-      m=iP*iT*iQ ! total number of layers
-
-      print*,'Temperature layers:  ',iT
-      print*,'Pressure layers:     ',iP
-      print*,'Mixing ratio layers: ',iQ
-      print*,'Total:               ',m
-
-      do mol=1,nmolec-1
-         print*,'Please enter vmr of ',molec_names(mol)
-         read(*,*) xconst(mol)
-      end do
-
-     
-      do i=1,iQ
-         do j=1,iP              
-            do k=1,iT
-
-               im = (i-1)*iP*iT + (j-1)*iT + k
-
-               alt(im)=0.D0
-               pres(im)=P(j)
-               temp(im)=T(k)
-
-               do mol=1,nmolec-1
-                  x(im,mol)=xconst(mol)*(1.0-Q(i))
-               enddo
-               x(im,nmolec)=Q(i)
-
-            enddo
-         enddo
-      enddo
-
-      open(12,file=composition_file(1:strlen(composition_file)))
-      string='Atmospheric composition input '
-      string=string(1:strlen(string))//' data file for planet:'
-      string=string(1:strlen(string))//' '//
-
-     &     planet_descriptor(1:strlen(planet_descriptor))
-      write(12,'(a)') string(1:strlen(string))
-      write(12,39) 'Number of atmospheric levels: ',m
-      write(12,34) 'Number of molecules: ',nmolec
-      write(12,*) 
-      string='      z (km)     /    P (atm)     /     T (K) *'
-
-
-      do mol=1,nmolec
-
-         molname=molec_names(mol)
-
-         nb=5
-         if (strlen(molname).eq.1) then
-            na=7
-         else if (strlen(molname).eq.2) then
-            na=6
-         else if (strlen(molname).eq.3) then
-            na=5
-         endif
-         s='*'
-         do ne=1,nb
-            s=s(1:strlen(s)-1)//spch//'*'
-         enddo
-         s=s(1:strlen(s)-1)//'/*'
-         do ne=1,na
-            s=s(1:strlen(s)-1)//spch//'*'
-         enddo
-         
-         string=string(1:strlen(string)-1)
-     &        //s(1:strlen(s)-1)
-     &        //'x['
-     &        //molname(1:strlen(molname))
-     &        //']*'
-      enddo
-      write(12,*) string(1:strlen(string)-1)
-
-      do im=1,m
-
-         do mol=1,nmolec
-            xlev(mol)=x(im,mol)
-         enddo                  ! mol
-
-         write(12,50) alt(im),pres(im),temp(im),(xlev(mol),mol=1,nmolec)
-
-      enddo
-
-      write(12,*) 
-      close(12)
-
-      write(*,*) 'Output file successfully generated:'
-      write(*,*) composition_file(1:strlen(composition_file))
-      
-      end
Index: trunk/LMDZ.GENERIC/utilities/max.inc
===================================================================
--- trunk/LMDZ.GENERIC/utilities/max.inc	(revision 1507)
+++ 	(revision )
@@ -1,23 +1,0 @@
-          integer Nline_mx,Nmol_max,Niso_max,Nbmx,Nkmx,Nmax,Nlpnb_mx
-	    integer Nxtab_mx,Nacc_mx,Nmax_tab,Nytab_mx,Nktab_mx
-	    integer nfiles_mx,Nproc_mx,Nchunk_mx,Nc_mx,Nmax_vt,Nppc_max
-	    integer Nti_mx
-          parameter(Nline_mx=2000000)
-          parameter(Nmol_max=42)
-          parameter(Niso_max=15)
-	  parameter(Nbmx=10000)
-	  parameter(Nkmx=100000)
-	  parameter(Nmax=2000)
-	  parameter(Nlpnb_mx=500000)
-	  parameter(Nxtab_mx=100000)
-	  parameter(Nacc_mx=1000)
-	  parameter(Nmax_tab=1000000)
-	    parameter(Nytab_mx=210)
-	    parameter(Nktab_mx=1000000)
-	    parameter(nfiles_mx=100)
-	    parameter(Nproc_mx=1000)
-	    parameter(Nchunk_mx=1000)
-	    parameter(Nc_mx=9)
-	    parameter(Nmax_vt=1000000)
-	    parameter(Nppc_max=10000)
-	    parameter(Nti_mx=100000)
Index: trunk/LMDZ.GENERIC/utilities/postkspectrum
===================================================================
--- trunk/LMDZ.GENERIC/utilities/postkspectrum	(revision 1507)
+++ 	(revision )
@@ -1,48 +1,0 @@
-#!/bin/bash
-
-DWORK_DIR=$RAD/CORRELATEDK/avec_robin_data
-KSPEC_DIR=$KSPECTRUM/kspec_run
-#BANDS_DIR=32x36
-BANDS_DIR=38x36
-
-### if hires_spectrum don't exist then create it, copy data from kspectrum
-DIRCHECK=hires_spectrum
-cd $DWORK_DIR/
-if [ -d $DIRCHECK ]; then
-    echo -en "Directory hires_spectrum already exists.\n"
-else
-    echo -en "Copying kspectrum data to its new home...\n"
-    cd $DWORK_DIR
-    mkdir hires_spectrum
-
-    mv $KSPEC_DIR/results/k* hires_spectrum/
-    cp $KSPEC_DIR/results/calculation_info.txt hires_spectrum/
-    cp $KSPEC_DIR/*in hires_spectrum/
-    cp $KSPEC_DIR/data/composition.in hires_spectrum/
-
-    # removing contents of 'optimizations' folder
-    rm -rf $KSPEC_DIR/optimizations
-    mkdir $KSPEC_DIR/optimizations
-
-    chmod 444 hires_spectrum/*
-fi
-
-cp ~/utilities/run_kmatrix .
-
-echo -en "\nNow you must submit the parallel job (>> llsubmit run_kmatrix)\n"
-
-exit 0
-
-#cp ~/utilities/generate_kmatrix.exe .
-### make kmatrix in the infrared and save it in a GCM-friendly format
-#cp $BANDS_DIR/narrowbands_IR.in narrowbands.in
-#generate_kmatrix.exe
-#mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_IR.dat 
-#rm narrowbands.in
-### make kmatrix in the visible and save it in a GCM-friendly format
-#cp $BANDS_DIR/narrowbands_VI.in narrowbands.in
-#generate_kmatrix.exe
-#mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_VI.dat 
-#rm narrowbands.in
-#rm generate_kmatrix.exe
-#exit 0
Index: trunk/LMDZ.GENERIC/utilities/prekspectrum
===================================================================
--- trunk/LMDZ.GENERIC/utilities/prekspectrum	(revision 1507)
+++ 	(revision )
@@ -1,28 +1,0 @@
-#!/bin/bash
-
-#DWORK_DIR=$RAD/CORRELATEDK/premier_essai
-
-DWORK_DIR=$PWD
-KSPEC_DIR=$KSPECTRUM/kspec_run
-
-cd $DWORK_DIR
-
-# test for presence of k* in kspec results directory
-FILECHECK=$KSPEC_DIR/results/k001
-if [ -f $FILECHECK ]; then
-    echo -en "There is already k-data in the kspec directory "
-    echo -en $KSPEC_DIR
-    echo -en ", delete or move it before starting!\n"
-    exit 1
-fi
-
-make_composition_gen.exe
-cp composition.in $KSPEC_DIR/data/
-
-echo -en "kspectrum working directory is:\n"
-echo -en $KSPEC_DIR
-echo -en "\nNow you must submit the parallel job (>> llsubmit run_kspec)\n"
-
-cd $KSPEC_DIR
-
-exit 0
Index: trunk/LMDZ.GENERIC/utilities/run_kmatrix
===================================================================
--- trunk/LMDZ.GENERIC/utilities/run_kmatrix	(revision 1507)
+++ 	(revision )
@@ -1,28 +1,0 @@
-# @ job_name=kmatrix
-# @ output = $(job_name).$(jobid)
-# @ error = $(output)
-# @ job_type = serial
-# @ class = B
-# @ queue
-
-BANDS_DIR=38x36
-
-cp $MOD/LMDZ.GENERIC/utilities/generate_kmatrix.exe .
-
-## make kmatrix in the infrared and save it in a GCM-friendly format
-cp $BANDS_DIR/narrowbands_IR.in narrowbands.in
-generate_kmatrix.exe > log.out
-mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_IR.dat 
-mv log.out $BANDS_DIR/logIR.out 
-rm narrowbands.in
-
-## make kmatrix in the visible and save it in a GCM-friendly format
-cp $BANDS_DIR/narrowbands_VI.in narrowbands.in
-generate_kmatrix.exe > log.out
-mv corrk_gcm.dat $BANDS_DIR/corrk_gcm_VI.dat 
-mv log.out $BANDS_DIR/logVI.out 
-rm narrowbands.in
-
-rm generate_kmatrix.exe
-
-exit 0
Index: trunk/LMDZ.GENERIC/utilities/startrestart
===================================================================
--- trunk/LMDZ.GENERIC/utilities/startrestart	(revision 1507)
+++ 	(revision )
@@ -1,6 +1,0 @@
-#!/bin/bash
-
-mv -f restart.nc start.nc
-mv -f restartfi.nc startfi.nc
-
-exit 0
Index: trunk/LMDZ.GENERIC/utilities/startrestart0
===================================================================
--- trunk/LMDZ.GENERIC/utilities/startrestart0	(revision 1507)
+++ 	(revision )
@@ -1,6 +1,0 @@
-#!/bin/bash
-
-mv -f restart.nc start0.nc
-mv -f restartfi.nc startfi0.nc
-
-exit 0
Index: trunk/LMDZ.GENERIC/utilities/strlen.for
===================================================================
--- trunk/LMDZ.GENERIC/utilities/strlen.for	(revision 1507)
+++ 	(revision )
@@ -1,10 +1,0 @@
-      integer function strlen(st)
-      integer i
-      character st*(*)
-      i=len(st)
-      do while (st(i:i).eq.' ')
-         i=i-1
-      enddo
-      strlen=i
-      return
-      end
