Index: trunk/LMDZ.TITAN/libf/phytitan/calc_cpp_mugaz.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calc_cpp_mugaz.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/calc_cpp_mugaz.F90	(revision 1648)
@@ -34,26 +34,15 @@
       do igas=1,ngasmx
 
-         if(igas.eq.vgas)then
-            ! ignore variable gas in cpp calculation
+         ! all values at 300 K from Engineering Toolbox
+         if(igas.eq.igas_N2)then
+            mugaz_c = mugaz_c + 28.01*gfrac(igas,nivref)
+         elseif(igas.eq.igas_H2)then
+            mugaz_c = mugaz_c + 2.01*gfrac(igas,nivref)
+         elseif(igas.eq.igas_CH4)then
+            mugaz_c = mugaz_c + 16.04*gfrac(igas,nivref)
          else
-            ! all values at 300 K from Engineering Toolbox
-            if(igas.eq.igas_N2)then
-               mugaz_c = mugaz_c + 28.01*gfrac(igas)
-            elseif(igas.eq.igas_H2)then
-               mugaz_c = mugaz_c + 2.01*gfrac(igas)
-            elseif(igas.eq.igas_CH4)then
-               mugaz_c = mugaz_c + 16.04*gfrac(igas)
-            elseif(igas.eq.igas_C2H6)then 
-               ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
-               mugaz_c = mugaz_c + 30.07*gfrac(igas)
-            elseif(igas.eq.igas_C2H2)then
-               ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
-               mugaz_c = mugaz_c + 26.04*gfrac(igas)
-            else
-               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
-               call abort
-            endif
+            print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
+            call abort
          endif
-
       enddo
 
@@ -61,33 +50,16 @@
       do igas=1,ngasmx
 
-         if(igas.eq.vgas)then
-            ! ignore variable gas in cpp calculation
+         ! all values at 300 K from Engineering Toolbox
+         if(igas.eq.igas_N2)then
+            cpp_c   = cpp_c   + 1.040*gfrac(igas,nivref)*28.01/mugaz_c
+         elseif(igas.eq.igas_H2)then
+            cpp_c   = cpp_c   + 14.31*gfrac(igas,nivref)*2.01/mugaz_c
+         elseif(igas.eq.igas_CH4)then
+            cpp_c   = cpp_c   + 2.226*gfrac(igas,nivref)*16.04/mugaz_c
          else
-            ! all values at 300 K from Engineering Toolbox
-            if(igas.eq.igas_N2)then
-               cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
-            elseif(igas.eq.igas_H2)then
-               cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
-            elseif(igas.eq.igas_CH4)then
-               cpp_c   = cpp_c   + 2.226*gfrac(igas)*16.04/mugaz_c
-            elseif(igas.eq.igas_C2H6)then
-               ! C2H6
-               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
-               cpp_c   = cpp_c   + 1.763*gfrac(igas)*30.07/mugaz_c
-            elseif(igas.eq.igas_C2H2)then
-               ! C2H2
-               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
-               cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
-            else
-               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
-               call abort
-            endif
+            print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
+            call abort
          endif
-
       enddo
-
-
-
-
 
       cpp_c = 1000.0*cpp_c
Index: trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90	(revision 1648)
@@ -26,5 +26,5 @@
 
       use radinc_h, only: L_NSPECTV
-      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
+      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep
       use gases_h
       use comcstfi_mod, only: g, mugaz
@@ -38,5 +38,5 @@
 
       logical typeknown
-      real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart
+      real*8 tauvar,tausum,tauwei,bwidth,bstart
       double precision df
 
@@ -52,8 +52,5 @@
       typeknown=.false.
       do igas=1,ngasmx
-         if(igas.eq.vgas)then
-            print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
-	 endif
-	 if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
+	 if(gfrac(igas,nivref).lt.5.e-2)then
             print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
             'as its mixing ratio is less than 0.05.' 
@@ -74,5 +71,5 @@
             endif
 
-            if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then
+            if(gfrac(igas,nivref).eq.1.0)then
                print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.'
                typeknown=.true.
@@ -92,6 +89,4 @@
          tausum = 0.0
          tauwei = 0.0
-         tausumvar = 0.0
-         tauweivar = 0.0
          bstart = 10000.0/BWNV(N+1)
          bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
@@ -100,7 +95,6 @@
 
             tauvar=0.0
-            tauvarvar=0.0
             do igas=1,ngasmx
-               if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
+               if(gfrac(igas,nivref).lt.5.e-2)then
                   ! ignore variable gas in Rayleigh calculation
                   ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
@@ -117,10 +111,5 @@
                endif
 
-               if(igas.eq.vgas) then
-                  tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas)
-                  tauvar=tauvar+0.0*0.0*gfrac(igas)
-               else
-                  tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas)
-               endif
+               tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas,nivref)
 
             enddo
@@ -130,10 +119,7 @@
             tauwei=tauwei+df
             tausum=tausum+tauvar*df
-            tauweivar=tauweivar+df
-            tausumvar=tausumvar+tauvarvar*df
          
          enddo
          TAURAY(N)=tausum/tauwei
-         TAURAYVAR(N)=tausumvar/tauweivar
          ! we multiply by scalep here (100) because plev, which is used in optcv,
          ! is in units of mBar, so we need to convert.
Index: trunk/LMDZ.TITAN/libf/phytitan/call_profilgases.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/call_profilgases.F90	(revision 1648)
+++ trunk/LMDZ.TITAN/libf/phytitan/call_profilgases.F90	(revision 1648)
@@ -0,0 +1,85 @@
+subroutine call_profilgases(nlayer)
+
+  use gases_h
+
+  implicit none
+
+
+  !============================================================================
+  !
+  !     Purpose
+  !     -------
+  !     Set atmospheric composition for Titan main gases (N2, CH4, H2) from :
+  !       - Interactive tracers if we run with chemistry
+  !       or
+  !       - CH4 vertical profile file (profile.def) at firstcall only
+  !
+  !     Author
+  !     ------
+  !     J. Vatant d'Ollone (2017)
+  !
+  !============================================================================
+
+
+  !--------------------------
+  ! 0. Declarations
+  ! -------------------------
+
+  integer, intent(in) :: nlayer
+  
+  
+  logical, save :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall)
+
+
+  integer ilay, ierr
+
+  
+  ! -----------------------------------------------
+  ! 1. First case : no interactive CH4 tracer
+  ! -----------------------------------------------
+
+  if (firstcall) then
+
+!$OMP MASTER
+     
+     ! Load CH4 vertical profile from file 'profile.def'
+     open(90,file='profile.def',status='old',form='formatted',iostat=ierr)
+
+     if (ierr.eq.0) then
+        write(*,*) "call_profilgases.F90: reading file profile.def"
+        read(90,*) ! header
+
+        write(*,*) "call_profilgases.F90: reading CH4 vertical profile..."
+
+        do ilay=1,nlayer
+           read(90,*,iostat=ierr) gfrac(igas_CH4,ilay)
+           if (ierr.ne.0) then
+              write(*,*) 'call_profilgases.F90: error reading CH4 vertical profile... aborting'
+              call abort
+           endif          
+        enddo        
+        
+        ! Then set H2 (fixed) and N2 (what remains)
+
+        gfrac(igas_H2,:) = 1.0E-3
+        gfrac(igas_N2,:) = 1.0 - ( gfrac(igas_CH4,:) + gfrac(igas_H2,:) )
+        
+     else
+        write(*,*) 'Cannot find required file "profile.def"'
+        call abort
+     endif
+
+     close(90)
+!$OMP END MASTER
+!$OMP BARRIER
+     
+     firstcall=.false.
+  endif ! if (firstcall)
+
+
+  ! -----------------------------------------------
+  ! 2. Second case : interactive CH4 tracer
+  ! -----------------------------------------------
+  
+end subroutine call_profilgases
Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1648)
@@ -1,5 +1,5 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
           albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    & 
-          tsurf,fract,dist_star,aerosol,muvar,                 &
+          tsurf,fract,dist_star,aerosol,                       &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
@@ -17,5 +17,5 @@
       USE tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: diurnal,tracer,nosurf,satval,        &
+      use callkeys_mod, only: diurnal,tracer,nosurf,        &
                               strictboundcorrk,specOLR
 
@@ -59,5 +59,4 @@
       REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
       REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
-      REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
       logical,intent(in) :: firstcall              ! Signals first call to physics.
       logical,intent(in) :: lastcall               ! Signals last call to physics.
@@ -133,5 +132,4 @@
       real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
       real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
-      real*8 qvar(L_LEVELS)                    ! Mixing ratio of variable component (mol/mol).
 
       ! Local aerosol optical properties for each column on RADIATIVE grid.
@@ -153,18 +151,7 @@
       character(len=10) :: tmp2
 
-      ! For fixed water vapour profiles.
-      integer i_var
-      real RH
-      real*8 pq_temp(nlayer)
-! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
-      real ptemp, Ttemp, qsat
-
       logical OLRz
       real*8 NFLUXGNDV_nu(L_NSPECTV)
-
-      ! Included by RW for runaway greenhouse 1D study.
-      real vtmp(nlayer)
-      REAL*8 muvarrad(L_LEVELS)
-      
+   
       ! Included by MT for albedo calculations.      
       REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
@@ -459,42 +446,4 @@
       else
          acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
-      endif
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!! Note by JL13 : In the following, some indices were changed in the interpolations,
-!!!                so that the model results are less dependent on the number of layers !
-!!!
-!!!           ---  The older versions are commented with the comment !JL13index  ---
-!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-     do k=1,L_LEVELS
-         qvar(k) = 1.0D-7
-     end do
-
-
-       ! IMPORTANT: Now convert from kg/kg to mol/mol.
-!       do k=1,L_LEVELS
-!          qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
-!       end do
-
-       DO l=1,nlayer
-          muvarrad(2*l)   = muvar(ig,nlayer+2-l)
-          muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
-       END DO
-      
-       muvarrad(1) = muvarrad(2)
-       muvarrad(2*nlayer+1)=muvar(ig,1)         
-      
-      ! Keep values inside limits for which we have radiative transfer coefficients !!!
-      if(L_REFVAR.gt.1)then ! (there was a bug here)
-         do k=1,L_LEVELS
-            if(qvar(k).lt.wrefvar(1))then
-               qvar(k)=wrefvar(1)+1.0e-8
-            elseif(qvar(k).gt.wrefvar(L_REFVAR))then
-               qvar(k)=wrefvar(L_REFVAR)-1.0e-8
-            endif
-         end do
       endif
 
@@ -574,16 +523,16 @@
             endif
          elseif(tlevrad(k).gt.tgasmax)then
-            print*,'Maximum temperature is outside the radiative'
-            print*,'transfer kmatrix bounds, exiting.'
-            print*,"k=",k," tlevrad(k)=",tlevrad(k)
-            print*,"tgasmax=",tgasmax
+!            print*,'Maximum temperature is outside the radiative'
+!            print*,'transfer kmatrix bounds, exiting.'
+!            print*,"k=",k," tlevrad(k)=",tlevrad(k)
+!            print*,"tgasmax=",tgasmax
             if (strictboundcorrk) then
               call abort
             else
-              print*,'***********************************************'
-              print*,'we allow model to continue with tlevrad=tgasmax'  
-              print*,'  ... we assume we know what you are doing ... '
-              print*,'  ... but do not let this happen too often ... '
-              print*,'***********************************************'
+!              print*,'***********************************************'
+!              print*,'we allow model to continue with tlevrad=tgasmax'  
+!              print*,'  ... we assume we know what you are doing ... '
+!              print*,'  ... but do not let this happen too often ... '
+!              print*,'***********************************************'
               !tlevrad(k)=tgasmax
             endif
@@ -607,16 +556,16 @@
             endif
          elseif(tmid(k).gt.tgasmax)then
-            print*,'Maximum temperature is outside the radiative'
-            print*,'transfer kmatrix bounds, exiting.'
-            print*,"k=",k," tmid(k)=",tmid(k)
-            print*,"tgasmax=",tgasmax
+!            print*,'Maximum temperature is outside the radiative'
+!            print*,'transfer kmatrix bounds, exiting.'
+!            print*,"k=",k," tmid(k)=",tmid(k)
+!            print*,"tgasmax=",tgasmax
             if (strictboundcorrk) then
               call abort
             else
-              print*,'***********************************************'
-              print*,'we allow model to continue with tmid=tgasmin'
-              print*,'  ... we assume we know what you are doing ... '
-              print*,'  ... but do not let this happen too often ... '
-              print*,'***********************************************'
+!              print*,'***********************************************'
+!              print*,'we allow model to continue with tmid=tgasmin'
+!              print*,'  ... we assume we know what you are doing ... '
+!              print*,'  ... but do not let this happen too often ... '
+!              print*,'***********************************************'
               tmid(k)=tgasmax
             endif
@@ -649,5 +598,5 @@
             call optcv(dtauv,tauv,taucumv,plevrad,                 &
                  qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
-                 tmid,pmid,taugsurf,qvar,muvarrad)
+                 tmid,pmid,taugsurf,gweight)
 
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
@@ -692,5 +641,5 @@
          call optci(plevrad,tlevrad,dtaui,taucumi,                  &
               qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
-              taugsurfi,qvar,muvarrad)
+              taugsurfi,gweight)
 
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
@@ -813,5 +762,4 @@
         IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
         IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
-        IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
         IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
 !$OMP END MASTER
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1648)
@@ -60,6 +60,5 @@
       real,save :: size_tropo
       real,save :: size_strato
-      real,save :: satval
-!$OMP THREADPRIVATE(size_tropo,size_strato,satval)
+!$OMP THREADPRIVATE(size_tropo,size_strato)
       real,save :: pceil
 !$OMP THREADPRIVATE(pceil)
Index: trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90	(revision 1648)
+++ trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90	(revision 1648)
@@ -0,0 +1,103 @@
+subroutine disr_haze(dz,press,wno,taeros,ssa,cbar)
+use datafile_mod, only: datadir
+IMPLICIT NONE
+
+real*8,intent(in)   :: dz,press,wno
+real*8,intent(inout):: taeros,ssa,cbar
+
+!---------------------------
+! Attention !!
+! taeros is the integrated extinction over the layer
+! (extinction * thickness of layer)
+!---------------------------
+
+integer             :: i,j,iw,ip,ierr
+real*8              :: wln,factw,factp
+integer,parameter   :: nbwl_PL=328,nblev_PL=162
+real*8,save         :: ext_PL(nblev_PL,nbwl_PL),ssa_PL(nblev_PL,nbwl_PL)
+real*8,save         :: asf_PL(nblev_PL,nbwl_PL)
+real*8,save         :: wl_PL(nbwl_PL),press_PL(nblev_PL)
+logical,save        :: firstcall=.true.
+character(len=15)   :: dummy
+
+if (firstcall) then
+   print*,"DISR HAZE - PANAYOTIS LAVVAS"
+
+
+! read PL table   
+! wl_PL in nm
+! press_PL in Pa
+   open(11,file=TRIM(datadir)//'/hazetable_PL_original.dat',status="old",iostat=ierr)
+   read(11,*) dummy,wl_PL
+   do i=1,nblev_PL
+     read(11,*) press_PL(i),ext_PL(i,:)  ! in cm-1
+   enddo
+   do i=1,nblev_PL
+     read(11,*) press_PL(i),ssa_PL(i,:)
+   enddo
+   do i=1,nblev_PL
+     read(11,*) press_PL(i),asf_PL(i,:)
+   enddo
+   close(11)
+   ! convert press_PL into millibar for comparison to press in the generic (wtf press in mbar but ok ... to be modified at some point anyway)
+   press_PL(:)=press_PL(:)*1E-2
+
+   firstcall=.false.
+endif
+
+! convert wno (in cm-1) into wln (nm)
+wln=1E7/wno
+
+! interpolate the needed values from the table
+
+iw=1
+do i=2,nbwl_PL
+  if(wln.gt.wl_PL(i)) then
+    iw=i
+  endif
+enddo
+
+ip=1
+do j=2,nblev_PL
+  if(press.lt.press_PL(j)) then
+    ip=j
+  endif
+enddo
+
+if(iw.ne.nbwl_PL) then
+  factw=(wln-wl_PL(iw))/(wl_PL(iw+1)-wl_PL(iw))
+endif
+if(ip.ne.nblev_PL) then
+  factp=(log10(press)-log10(press_PL(ip)))/(log10(press_PL(ip+1))-log10(press_PL(ip)))
+else
+!  factp=(log10(press)-log10(press_PL(nblev_PL)))/(log10(1.e-10)-log10(press_PL(nblev_PL)))
+  factp=0.
+endif
+
+if ((iw.ne.nbwl_PL).and.(ip.ne.nblev_PL)) then
+  taeros= ext_PL(ip,iw)  *(1-factw)*(1-factp)+ext_PL(ip+1,iw)  *(1-factw)*factp &
+         +ext_PL(ip,iw+1)*factw    *(1-factp)+ext_PL(ip+1,iw+1)*factw    *factp
+  ssa   = ssa_PL(ip,iw)  *(1-factw)*(1-factp)+ssa_PL(ip+1,iw)  *(1-factw)*factp &
+         +ssa_PL(ip,iw+1)*factw    *(1-factp)+ssa_PL(ip+1,iw+1)*factw    *factp
+  cbar  = asf_PL(ip,iw)  *(1-factw)*(1-factp)+asf_PL(ip+1,iw)  *(1-factw)*factp &
+         +asf_PL(ip,iw+1)*factw    *(1-factp)+asf_PL(ip+1,iw+1)*factw    *factp
+else if ((iw.eq.nbwl_PL).and.(ip.ne.nblev_PL)) then
+  taeros= ext_PL(ip,iw)*(1-factp)+ext_PL(ip+1,iw)*factp 
+  ssa   = ssa_PL(ip,iw)*(1-factp)+ssa_PL(ip+1,iw)*factp 
+  cbar  = asf_PL(ip,iw)*(1-factp)+asf_PL(ip+1,iw)*factp
+else if ((iw.ne.nbwl_PL).and.(ip.eq.nblev_PL)) then
+  taeros= ext_PL(ip,iw)  *(1-factw)*(1-factp) &
+         +ext_PL(ip,iw+1)*factw    *(1-factp)
+  ssa   = ssa_PL(ip,iw)  *(1-factw)*(1-factp) &
+         +ssa_PL(ip,iw+1)*factw    *(1-factp)
+  cbar  = asf_PL(ip,iw)  *(1-factw)*(1-factp) &
+         +asf_PL(ip,iw+1)*factw    *(1-factp)
+else if ((iw.eq.nbwl_PL).and.(ip.eq.nblev_PL)) then
+  taeros= ext_PL(ip,iw)*(1-factp) 
+  ssa   = ssa_PL(ip,iw)*(1-factp) 
+  cbar  = asf_PL(ip,iw)*(1-factp)
+endif
+
+taeros=taeros*dz*1.E2 ! ext in cm-1 * thickness in m * 1E2
+
+end subroutine disr_haze
Index: trunk/LMDZ.TITAN/libf/phytitan/gases_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/gases_h.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/gases_h.F90	(revision 1648)
@@ -3,26 +3,29 @@
       implicit none
 
-!======================================================================C
+!============================================================================================C
 !
 !     gases_h    
 !
-!     A F90-allocatable version for gases.h -- AS 12/2011
+!     * A F90-allocatable version for gases.h -- AS 12/2011
 !
-!======================================================================C
+!     * Titan's version : J. Vatant d'Ollone (2017)
+!                  * gfrac is now 2D (altitude-dependent for the CIA-calculations) 
+!                  * Reference level (nivref) is added for the cpp_mu and rayleigh routines
+!===========================================================================================C
 
       ! Set and allocated in su_gases.F90
       integer :: ngasmx
-      integer :: vgas
-      character*20,allocatable,DIMENSION(:) :: gnom ! name of the gas, read by master
-      real,allocatable,DIMENSION(:) :: gfrac
+      integer :: nivref
+      character*20,allocatable,DIMENSION(:) :: gnom ! name of the gas, set by master
+      real,allocatable,DIMENSION(:,:) :: gfrac
 
       ! in analogy with tracer.h ...
-      integer :: igas_H2
+      
       integer :: igas_N2
       integer :: igas_CH4
-      integer :: igas_C2H2
-      integer :: igas_C2H6
-!!$OMP THREADPRIVATE(ngasmx,vgas,gnom,gfrac,&
-!	!$OMP igas_H2,igas_N2,igas_CH4,igas_C2H2,igas_C2H6)
+      integer :: igas_H2
+      
+!!$OMP THREADPRIVATE(ngasmx,nivref,gnom,gfrac,&
+!	!$OMP igas_H2,igas_N2,igas_CH4)
 
       end module gases_h
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1648)
@@ -439,9 +439,4 @@
 
 !=================================
-
-     write(*,*)"What is the saturation % of the variable species?"
-     satval=0.8
-     call getin_p("satval",satval)
-     write(*,*)" satval = ",satval
 
      write(*,*) "Gravitationnal sedimentation ?"
@@ -477,7 +472,10 @@
 !           mugaz=8.314*1000./pr
      endif ! of if (force_cpp)
-     call su_gases
+     
+     
+     call su_gases(nlayer,tracer)     
      call calc_cpp_mugaz
-
+     
+     
      PRINT*,'--------------------------------------------'
      PRINT*
@@ -519,5 +517,4 @@
   ENDDO
 
-  ! initialize variables in radinc_h
   call ini_radinc_h(nlayer)
  
Index: trunk/LMDZ.TITAN/libf/phytitan/interpolateCH4CH4.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/interpolateCH4CH4.F90	(revision 1648)
+++ trunk/LMDZ.TITAN/libf/phytitan/interpolateCH4CH4.F90	(revision 1648)
@@ -0,0 +1,128 @@
+subroutine interpolateCH4CH4(wn,temp,pres,abcoef,firstcall,ind)
+
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Calculates the CH4-CH4 CIA opacity, using a lookup table from
+  !     HITRAN (2011)
+  !
+  !     Authors
+  !     -------
+  !     J. Vatant (2016) based on 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)
+  integer :: ind
+
+  ! output
+  double precision abcoef             ! absorption coefficient (m^-1)
+
+  integer nS,nT
+  parameter(nS=1018)
+  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 !read by master
+
+  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 CH4-CH4 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 CH4-CH4 continuum from HITRAN database...'
+
+     !     1.1 Open the ASCII files
+     dt_file=TRIM(datadir)//'/continuum_data/CH4-CH4_2011.cia'
+
+!$OMP MASTER
+     open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+     if (ios.ne.0) then        ! file not found
+        write(*,*) 'Error from interpolateCH4CH4'
+        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/CH4-CH4_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 interpolateCH4CH4.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)
+!$OMP END MASTER
+!$OMP BARRIER
+
+     print*,'interpolateCH4CH4: At wavenumber ',wn,' cm^-1'
+     print*,'   temperature ',temp,' K'
+     print*,'   pressure ',pres,' Pa'
+
+  endif
+     call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
+
+!     print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
+!     print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
+
+     abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
+!     abcoef=0.
+
+!     print*,'We have ',amagat,' amagats of CH4'
+!     print*,'So the absorption is ',abcoef,' m^-1'
+
+!     unlike for Rayleigh scattering, we do not currently weight by the BB function
+!     however our bands are normally thin, so this is no big deal.
+
+
+  return
+end subroutine interpolateCH4CH4
+
Index: trunk/LMDZ.TITAN/libf/phytitan/interpolateN2CH4.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/interpolateN2CH4.F90	(revision 1648)
+++ trunk/LMDZ.TITAN/libf/phytitan/interpolateN2CH4.F90	(revision 1648)
@@ -0,0 +1,133 @@
+subroutine interpolateN2CH4(wn,temp,presN2,presCH4,abcoef,firstcall,ind)
+
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Calculates the N2-CH4 CIA opacity, using a lookup table from
+  !     HITRAN (2011)
+  !
+  !     Authors
+  !     -------
+  !     J. Vatant (2016) based on 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 presCH4            ! CH4 partial pressure   (Pascals)
+
+  ! output
+  double precision abcoef             ! absorption coefficient (m^-1)
+
+  integer nS,nT
+  parameter(nS=1407)
+  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 amagatCH4
+  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 !read by master
+
+  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
+  integer ind
+
+  if(temp.gt.400)then
+     print*,'Your temperatures are too high for this N2-CH4 CIA dataset.'
+     print*,'Please run mixed N2-CH4 atmospheres below T = 400 K.'      
+     stop
+  endif
+
+  amagatN2 = (273.15/temp)*(presN2/101325.0)
+  amagatCH4 = (273.15/temp)*(presCH4/101325.0)
+
+  if(firstcall)then ! called by sugas_corrk only
+     print*,'----------------------------------------------------'
+     print*,'Initialising N2-CH4 continuum from HITRAN database...'
+
+     !     1.1 Open the ASCII files
+     dt_file=TRIM(datadir)//'/continuum_data/N2-CH4_2011.cia'
+
+!$OMP MASTER
+     open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+     if (ios.ne.0) then        ! file not found
+        write(*,*) 'Error from interpolateN2CH4'
+        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-CH4_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 interpolateN2CH4.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)
+!$OMP END MASTER
+!$OMP BARRIER
+
+     print*,'interpolateN2CH4: At wavenumber ',wn,' cm^-1'
+     print*,'   temperature                  ',temp,' K'
+     print*,'   N2 partial pressure          ',presN2,' Pa'
+     print*,'   and CH4 partial pressure     ',presCH4,' Pa'
+
+  endif
+
+  call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
+
+!     print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
+!     print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
+
+  abcoef=abcoef*losch**2*100.0*amagatN2*amagatCH4 ! convert to m^-1
+
+!     print*,'We have ',amagatN2,' amagats of N2'
+!     print*,'and     ',amagatCH4,' amagats of CH4'
+!     print*,'So the absorption is ',abcoef,' m^-1'
+
+
+!  unlike for Rayleigh scattering, we do not currently weight by the BB function
+!  however our bands are normally thin, so this is no big deal.
+
+
+  return
+end subroutine interpolateN2CH4
+
Index: trunk/LMDZ.TITAN/libf/phytitan/interpolateN2N2.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/interpolateN2N2.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/interpolateN2N2.F90	(revision 1648)
@@ -78,5 +78,5 @@
         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.'
+        write(*,*) 'Also check that the continuum data continuum_data/N2-N2_2011.cia is there.'
         call abort
      else
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1648)
@@ -1,10 +1,10 @@
 subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
      QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
-     TMID,PMID,TAUGSURF,QVAR,MUVAR)
+     TMID,PMID,TAUGSURF,GWEIGHT)
 
   use radinc_h
-  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
+  use radcommon_h, only: gasi,tlimit,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
   use gases_h
-  use comcstfi_mod, only: g, r, mugaz
+  use comcstfi_mod, only: g, r
   use callkeys_mod, only: continuum,graybody
   implicit none
@@ -49,4 +49,17 @@
   real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
 
+  ! Titan customisation
+  ! J. Vatant d'Ollone (2016)
+  real*8 GWEIGHT(L_NGAUSS)
+  real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI)
+  real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI)
+  real*8 SSA_T(L_LEVELS+1,L_NSPECTI)
+  real*8 ASF_T(L_LEVELS+1,L_NSPECTI)
+  real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
+  real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
+
+  CHARACTER*2  str2
+  ! ==========================
+
   integer L, NW, NG, K, LK, IAER
   integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
@@ -60,9 +73,6 @@
   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
@@ -73,5 +83,5 @@
   !real*8 rho !! see test below
 
-  integer igas, jgas
+  integer igas, jgas, ilay
 
   integer interm
@@ -95,10 +105,8 @@
 
      ! if we have continuum opacities, we need dz
-      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
-      U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
-	    !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
-
-     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
-          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
+      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))
+      U(k)  = Cmk*DPR(k)     ! only Cmk line in optci.F
+      
+     call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
 
      do LK=1,4
@@ -120,4 +128,6 @@
      do K=2,L_LEVELS
 
+        ilay = k / 2 ! int. arithmetic => gives the gcm layer index
+        
 ! continuum absorption
         DCONT = 0.0d0 
@@ -129,9 +139,5 @@
            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
+              p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
 
               dtemp=0.0d0
@@ -151,9 +157,28 @@
                  ! then cross-interactions with other gases
                  do jgas=1,ngasmx
-                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
                     dtempc  = 0.0d0
                     if(jgas.eq.igas_N2)then 
                        interm = indi(nw,igas,jgas)
                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indi(nw,igas,jgas) = interm
+                    endif
+                    dtemp = dtemp + dtempc
+                 enddo
+
+               elseif(igas.eq.igas_CH4)then
+
+                 ! first do self-induced absorption
+                 interm = indi(nw,igas,igas)
+                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
+                 indi(nw,igas,igas) = interm
+
+                 ! then cross-interactions with other gases
+                 do jgas=1,ngasmx
+                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
+                    dtempc  = 0.0d0
+                    if(jgas.eq.igas_N2)then 
+                       interm = indi(nw,igas,jgas)
+                       call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
                        indi(nw,igas,jgas) = interm
                     endif
@@ -185,36 +210,17 @@
 	DAERO=SUM(TAEROS(K,NW,1:naerkind))
 
+!================= Titan customisation ========================================
+        call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+! =============================================================================
+
+
         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
+           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)
 
            ! Interpolate the gaseous k-coefficients to the requested T,P values
@@ -228,5 +234,6 @@
            DTAUKI(K,nw,ng) = TAUGAS    & 
                              + DCONT   & ! For parameterized continuum absorption
-			     + DAERO     ! For aerosol absorption
+			     + DAERO   & ! For aerosol absorption
+			     + DHAZE_T(K,NW)  ! For Titan haze
 
         end do
@@ -238,5 +245,6 @@
         DTAUKI(K,nw,ng) = 0.d0      & 
                           + DCONT   & ! For parameterized continuum absorption
-	                  + DAERO     ! For aerosol absorption
+	                  + DAERO   & ! For aerosol absorption
+	                  + DHAZE_T(K,NW)     ! For Titan Haze
 
      end do
@@ -244,7 +252,10 @@
 
   DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0
+  
+  
   !=======================================================================
   !     Now the full treatment for the layers, where besides the opacity
   !     we need to calculate the scattering albedo and asymmetry factors
+  ! ======================================================================
 
   do iaer=1,naerkind
@@ -256,17 +267,34 @@
   end do
   
+  ! Haze scattering
   DO NW=1,L_NSPECTI
-     DO L=1,L_NLAYRAD
+    DO K=2,L_LEVELS+1
+      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
+    ENDDO
+  ENDDO
+
+  DO NW=1,L_NSPECTI
+     DO L=1,L_NLAYRAD-1
         K              = 2*L+1
-        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
+        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
+                      + DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
      END DO ! L vertical loop
+     
+     !last level
+     L              = L_NLAYRAD
+     K              = 2*L+1
+     
+     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + DHAZES_T(K,NW)
+     
   END DO                    ! NW spectral loop
   
+! ======================================================================
 
   DO NW=1,L_NSPECTI
      NG = L_NGAUSS
-     DO L=1,L_NLAYRAD
+     DO L=1,L_NLAYRAD-1
 
         K              = 2*L+1
+
         DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
 
@@ -278,4 +306,5 @@
                    GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
            end do
+           atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
            WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
         else
@@ -291,18 +320,42 @@
 
      END DO ! L vertical loop
+     
+     !     No vertical averaging on bottom layer
+
+     L = L_NLAYRAD
+     K = 2*L+1
+     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
+     
+        atemp = 0.
+        if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
+           do iaer=1,naerkind
+              atemp = atemp +                                     &
+                   GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
+           end do
+           atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
+           WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
+        else
+           WBARI(L,nw,ng) = 0.0D0
+           DTAUI(L,NW,NG) = 1.0D-9
+        endif
+
+        if(btemp(L,nw) .GT. 0.0d0) then
+           cosbi(L,NW,NG) = atemp/btemp(L,nw)
+        else
+           cosbi(L,NW,NG) = 0.0D0
+        end if
 
      ! 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
+           DO L=1,L_NLAYRAD-1
               K              = 2*L+1
               DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
 
               if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
-
-                 WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
-
+                WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
               else
                  WBARI(L,nw,ng) = 0.0D0
@@ -312,4 +365,18 @@
               cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
            END DO ! L vertical loop
+           
+            !     No vertical averaging on bottom layer
+
+            L = L_NLAYRAD
+            K = 2*L+1
+            DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
+            if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
+              WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
+            else
+               WBARI(L,nw,ng) = 0.0D0
+               DTAUI(L,NW,NG) = 1.0D-9
+            endif
+            cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
+           
         END IF
 
@@ -334,16 +401,22 @@
   ! 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)
-  
-!  print*,'WBARI'
-!  print*,WBARI
-!  print*,'DTAUI'
-!  print*,DTAUI
-!  call abort
-  
+
+!  Titan's outputs (J.V.O, 2016)===============================================
+!      do l=1,L_NLAYRAD
+!         do nw=1,L_NSPECTI
+!          INT_DTAU(L,NW) = 0.0d+0
+!            DO NG=1,L_NGAUSS
+!               INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtaui(L,nw,ng)*gweight(NG)
+!            enddo
+!         enddo
+!      enddo
+
+!	do nw=1,L_NSPECTI
+!          write(str2,'(i2.2)') nw
+!	  call writediagfi(1,'kgi'//str2,'Gaz extinction coefficient IR band '//str2,'m-1',1,int_dtau(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))
+!          call writediagfi(1,'khi'//str2,'Haze extinction coefficient IR band '//str2,'m-1',1,k_haze(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))	
+!	enddo  
+  
+! ==============================================================================  
 
   return
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1648)
@@ -1,10 +1,10 @@
 SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
      QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
-     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
+     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,GWEIGHT)
 
   use radinc_h
-  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
+  use radcommon_h, only: gasv, tlimit, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
   use gases_h
-  use comcstfi_mod, only: g, r, mugaz
+  use comcstfi_mod, only: g, r
   use callkeys_mod, only: continuum,graybody,callgasvis
 
@@ -55,4 +55,17 @@
   real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
 
+  ! Titan customisation
+  ! J. Vatant d'Ollone (2016)
+  real*8 GWEIGHT(L_NGAUSS)
+  real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI)
+  real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI)
+  real*8 SSA_T(L_LEVELS+1,L_NSPECTI)
+  real*8 ASF_T(L_LEVELS+1,L_NSPECTI)
+  real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
+  real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
+  
+  CHARACTER*2  str2
+  ! ==========================
+
   integer L, NW, NG, K, LK, IAER
   integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
@@ -69,8 +82,5 @@
   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
@@ -82,5 +92,5 @@
   real*8 dz(L_LEVELS)
 
-  integer igas, jgas
+  integer igas, jgas, ilay
 
   integer interm
@@ -107,11 +117,8 @@
      ! if we have continuum opacities, we need dz
 
-      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
-      U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
-	  !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
-     
-
-     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
-          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
+      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))
+      U(k)  = Cmk*DPR(k)     ! only Cmk line in optcv.F     
+
+     call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
 
      do LK=1,4
@@ -133,7 +140,9 @@
      end do                    ! levels
   end do
-  
+
   !     we ignore K=1...
   do K=2,L_LEVELS
+  
+     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
 
      do NW=1,L_NSPECTV
@@ -153,9 +162,5 @@
            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
+              p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
 
               dtemp=0.0
@@ -176,5 +181,5 @@
                  ! then cross-interactions with other gases
                  do jgas=1,ngasmx
-                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
                     dtempc  = 0.0
                     if(jgas.eq.igas_N2)then 
@@ -187,4 +192,23 @@
                  enddo
 
+               elseif(igas.eq.igas_CH4)then
+
+                 ! first do self-induced absorption
+                 interm = indv(nw,igas,igas)
+                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
+                 indv(nw,igas,igas) = interm
+
+                 ! then cross-interactions with other gases
+                 do jgas=1,ngasmx
+                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
+                    dtempc  = 0.0
+                    if(jgas.eq.igas_N2)then 
+                       interm = indv(nw,igas,jgas)
+                       call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indv(nw,igas,jgas) = interm
+                    endif
+                    dtemp = dtemp + dtempc
+                 enddo
+
               endif
 
@@ -197,36 +221,16 @@
         endif
 
+!================= Titan customisation ========================================
+ 	call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+! =============================================================================
+
         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
+           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)
 
            ! Interpolate the gaseous k-coefficients to the requested T,P values
@@ -240,5 +244,6 @@
            DTAUKV(K,nw,ng) = TAUGAS & 
                              + TRAYAER & ! TRAYAER includes all scattering contributions
-                             + DCONT ! For parameterized continuum aborption
+                             + DCONT    & ! For parameterized continuum aborption
+			     + DHAZE_T(K,NW)  ! For Titan haze
 
         end do
@@ -248,5 +253,6 @@
 
         NG              = L_NGAUSS
-        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
+        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT  & ! For parameterized continuum absorption
+			      + DHAZE_T(K,NW)  ! For Titan haze
 
         do iaer=1,naerkind
@@ -261,5 +267,6 @@
   !     Now the full treatment for the layers, where besides the opacity
   !     we need to calculate the scattering albedo and asymmetry factors
-
+  ! ======================================================================
+  
   do iaer=1,naerkind
     DO NW=1,L_NSPECTV
@@ -269,10 +276,24 @@
     ENDDO
   end do
+  
+  ! Haze scattering
+  DO NW=1,L_NSPECTV
+    DO K=2,L_LEVELS+1
+      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
+    ENDDO
+  ENDDO
+
 
   DO NW=1,L_NSPECTV
      DO L=1,L_NLAYRAD-1
         K              = 2*L+1
-	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind))
-        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
+        
+	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) &
+                    + SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind)) &
+                    + ASF_T(K,NW)*DHAZES_T(K,NW)
+                    
+        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
+                    + DHAZES_T(K,NW)
+        
 	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
 	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
@@ -283,6 +304,9 @@
      L              = L_NLAYRAD
      K              = 2*L+1
-     atemp(L,NW)    = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
-     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
+     
+     atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) &
+                 + ASF_T(K,NW)*DHAZES_T(K,NW)
+     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) &
+                 + DHAZES_T(K,NW)
      ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW)
      btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
@@ -292,4 +316,6 @@
   END DO                    ! NW spectral loop
 
+! ===========================================================================================
+
   DO NG=1,L_NGAUSS
     DO NW=1,L_NSPECTV
@@ -309,4 +335,5 @@
 
         WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
+
      END DO                 ! NW spectral loop
   END DO                    ! NG Gauss loop
@@ -329,4 +356,23 @@
 
 
+!  Titan's outputs (J.V.O, 2016)===============================================
+!      do l=1,L_NLAYRAD
+!         do nw=1,L_NSPECTV
+!          INT_DTAU(L,NW) = 0.0d+0
+!            DO NG=1,L_NGAUSS
+!               INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtauv(L,nw,ng)*gweight(NG) 
+!            enddo
+!         enddo
+!      enddo
+
+!	do nw=1,L_NSPECTV
+!          write(str2,'(i2.2)') nw
+!	  call writediagfi(1,'kgv'//str2,'Gaz extinction coefficient VI band '//str2,'m-1',1,int_dtau(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))
+!          call writediagfi(1,'khv'//str2,'Haze extinction coefficient VI band '//str2,'m-1',1,k_haze(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))	
+!	enddo  
+
+! ==============================================================================  
+
+
   return
 
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1648)
@@ -15,5 +15,4 @@
  
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
-      use gases_h, only: gnom, gfrac
       use radcommon_h, only: sigma, glat, grav, BWNV
       use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
@@ -23,5 +22,5 @@
       use geometry_mod, only: latitude, longitude, cell_area
       USE comgeomfi_h, only: totarea, totarea_planet
-      USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
+      USE tracer_h, only: noms, mmol, radius, qext, &
                           alpha_lift, alpha_devil, qextrhor, &
                           igcm_dustbin
@@ -353,6 +352,4 @@
 !$OMP THREADPRIVATE(qsurf_hist)
 
-      real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
-
       real,allocatable,dimension(:,:,:),save :: reffrad  ! Aerosol effective radius (m). By RW
       real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW
@@ -660,11 +657,11 @@
 ! II.a Call correlated-k radiative transfer scheme
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- 
-               muvar(1:ngrid,1:nlayer+1)=mugaz
+
+               call call_profilgases(nlayer)
 
                ! standard callcorrk
                call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
-                              tsurf,fract,dist_star,aerosol,muvar,                &
+                              tsurf,fract,dist_star,aerosol,                      &
                               zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
                               fluxsurfabs_sw,fluxtop_lw,                          &
Index: trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 1648)
@@ -7,5 +7,5 @@
 !
 !                             radcommon.h
-!
+!v
 !----------------------------------------------------------------------C
 !
@@ -36,6 +36,4 @@
 !     TAURAY     - Array (NSPECTV elements) of the pressure-independent
 !                  part of Rayleigh scattering optical depth.
-!     TAURAYVAR  - Array (NSPECTV elements) of the pressure-independent
-!                  part of Rayleigh scattering optical depth for the variable gas.
 !     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
 !                  each temperature, pressure, and spectral interval
@@ -63,8 +61,8 @@
       REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi
       REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in setspv
-      REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV)
+      REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV)
 !$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,&
 	!$OMP WNOV,DWNV,WAVEV,&
-	!$OMP STELLARF,TAURAY,TAURAYVAR)
+	!$OMP STELLARF,TAURAY)
 
       REAL*8 blami(L_NSPECTI+1)
@@ -79,10 +77,10 @@
       !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011  
       REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv
-      REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, WREFVAR, PFGASREF
+      REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, PFGASREF
       real*8 FZEROI(L_NSPECTI)
       real*8 FZEROV(L_NSPECTV)
       real*8 pgasmin, pgasmax
       real*8 tgasmin, tgasmax
-!$OMP THREADPRIVATE(gasi,gasv,&  !wrefvar,pgasref,tgasref,pfgasref read by master in sugas_corrk
+!$OMP THREADPRIVATE(gasi,gasv,&  !pgasref,tgasref,pfgasref read by master in sugas_corrk
 	!$OMP FZEROI,FZEROV)     !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
 
Index: trunk/LMDZ.TITAN/libf/phytitan/su_gases.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/su_gases.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/su_gases.F90	(revision 1648)
@@ -1,3 +1,3 @@
-subroutine su_gases
+subroutine su_gases(nlayer,tracer)
 
   use gases_h
@@ -5,6 +5,9 @@
   implicit none
 
-  integer igas, ierr, count
-
+  integer,intent(in) :: nlayer
+  logical,intent(in) :: tracer
+  
+  integer igas, ierr
+  
   !==================================================================
   !     
@@ -17,85 +20,67 @@
   !     R. Wordsworth (2011)
   !     Allocatable arrays by A. Spiga (2011)
-  !     
+  !     Titan's version : J. Vatant d'Ollone (2017)
+  !              * force igas_X labels and def gnom for N2,CH4,H2
+  !              * get rid of variable species ( enrichment dimension will be set in corrk routines )
+  !
   !==================================================================
 
 !$OMP MASTER
-  ! load gas names from file 'gases.def'
+
+
+  ngasmx   = 3 ! N2, CH4, H2
+
+
+  ! load reference level and reference molar fractions 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
+     read(90,*) ! header
+
+     ! We allocate gfrac and we set gas molar fractions for the reference level only.
+     ! This will be useful for the cpp_mu and rayleigh routines
+     ! Other gas molar fractions are now set in routine callprofilgases
+     
+     write(*,*) 'sugases.F90: allocating and reading gas molar fractions from reference level in gases.def...'
+     
+     if(.not.allocated(gfrac)) allocate(gfrac(ngasmx,nlayer))
+     
+     read(90,*,iostat=ierr) nivref     
      if (ierr.ne.0) then
-        write(*,*) "sugases.F90: error reading number of gases"
+        write(*,*) "sugases.F90: error reading reference level"
         write(*,*) "   (first line of gases.def) "
         call abort
      endif
-
-     print*,ngasmx, " gases found in gases.def. Allocating names and molar fractions..."
-
-     if (.not.allocated(gnom)) allocate(gnom(ngasmx))
+     
+     print*, "layer", nivref, "is reference level found in gases.def ..."
+     
      do igas=1,ngasmx
-        read(90,*,iostat=ierr) gnom(igas)
+        read(90,*,iostat=ierr) gfrac(igas,nivref)
         if (ierr.ne.0) then
-           write(*,*) 'sugases.F90: error reading gas names in gases.def...'
+           write(*,*) 'sugases.F90: error reading reference gas molar fractions in gases.def... aborting'
            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
 
-        ! 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
+     ! We force gnom = (N2, CH4, H2) and igas_X for Titan
+     
+     igas_N2  = 1
+     igas_CH4 = 2
+     igas_H2  = 3
 
 
-     ! assign the 'igas_X' labels
-     count=0
-     do igas=1,ngasmx
-        if (trim(gnom(igas)).eq."H2_" .or. trim(gnom(igas)).eq."H2") then
-           igas_H2=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."N2_" .or. trim(gnom(igas)).eq."N2") then
-           igas_N2=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."CH4") then
-           igas_CH4=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."C2H6") then
-           igas_C2H6=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."C2H2") then
-           igas_C2H2=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
-
+     if (.not.allocated(gnom)) allocate(gnom(ngasmx))
+     gnom(igas_N2)  = "N2_"
+     gnom(igas_CH4) = "CH4"
+     gnom(igas_H2)  = "H2_"
+     
+     
   else
      write(*,*) 'Cannot find required file "gases.def"'
      call abort
   endif
+  
   close(90)
 !$OMP END MASTER
Index: trunk/LMDZ.TITAN/libf/phytitan/sugas_corrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/sugas_corrk.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/sugas_corrk.F90	(revision 1648)
@@ -24,5 +24,5 @@
       use radcommon_h, only : tgasref,tgasmin,tgasmax
       use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
-      use radcommon_h, only : wrefvar,WNOI,WNOV
+      use radcommon_h, only : WNOI,WNOV
       use datafile_mod, only: datadir
       use comcstfi_mod, only: mugaz
@@ -44,5 +44,4 @@
       ! ALLOCATABLE ARRAYS -- AS 12/2011
       REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8 	!read by master
-      character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for check with gnom, read by master
 
       real*8 x, xi(4), yi(4), ans, p
@@ -79,10 +78,4 @@
       read(111,*) ngas
 
-      if(ngas.ne.ngasmx)then
-         print*,'Number of gases in radiative transfer data (',ngas,') does not', &
-                'match that in gases.def (',ngasmx,'), exiting.'
-         call abort
-      endif 
-
       if(ngas.gt.5 .or. ngas.lt.1)then
          print*,ngas,' species in database [',               &
@@ -90,40 +83,7 @@
                 '], radiative code cannot handle this.'
          call abort
-      endif 
-
-      ! dynamically allocate gastype and read from Q.dat
-      IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) )
-
-      do igas=1,ngas
-         read(111,*) gastype(igas)
-         print*,'Gas ',igas,' is ',gastype(igas)
-      enddo
-
-      ! get array size, load the coefficients
-      open(111,file=TRIM(file_path),form='formatted')
-      read(111,*) L_REFVAR
-      IF( .NOT. ALLOCATED( wrefvar ) ) ALLOCATE( WREFVAR(L_REFVAR) )
-      read(111,*) wrefvar
-      close(111)
-
-      ! Check that gastype and gnom match
-      do igas=1,ngas
-         print*,'Gas ',igas,' is ',trim(gnom(igas))
-         if (trim(gnom(igas)).ne.trim(gastype(igas))) then
-            print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', &
-                 'match that in gases.def (',trim(gnom(igas)),'), exiting. You should compare ', &
-                 'gases.def with Q.dat in your radiative transfer directory.' 
-            call abort
-         endif
-      enddo
-      print*,'Confirmed gas match in radiative transfer and gases.def!'
-
-      ! display the values
-      print*,'Variable gas volume mixing ratios:'
-      do n=1,L_REFVAR
-         !print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
-         print*,n,'.',wrefvar(n),' mol/mol'
-      end do
-      print*,''
+      endif
+      
+      L_REFVAR = 1 ! JVO 2017 : set to 1 to keep the code running until the new variable species treatment
 
 !=======================================================================
@@ -628,6 +588,19 @@
             enddo
 
-         endif
-
+         elseif (igas .eq. igas_CH4) then
+
+            ! first do self-induced absorption
+            dummy = -9999
+            call interpolateCH4CH4(200.D+0,200.D+0,7500.D+0,testcont,.true.,dummy)
+            ! then cross-interactions with other gases
+            do jgas=1,ngasmx
+               if (jgas .eq. igas_N2) then
+                  dummy = -9999
+                  call interpolateN2CH4(200.D+0,250.0D+0,100000.D+0,5000.D+0,testcont,.true.,dummy)
+               endif
+            enddo
+
+         endif  
+         
       enddo
       endif
@@ -645,5 +618,4 @@
       IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
       IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
-      IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
 !$OMP END MASTER
 !$OMP BARRIER
Index: trunk/LMDZ.TITAN/libf/phytitan/tpindex.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tpindex.F	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/tpindex.F	(revision 1648)
@@ -1,4 +1,3 @@
-      subroutine tpindex(pw,tw,qvar,pref,tref,wrefvar,LCOEF,MT,MP,
-     &     NVAR,wratio)
+      subroutine tpindex(pw,tw,pref,tref,LCOEF,MT,MP)
 
 !==================================================================
@@ -6,5 +5,5 @@
 !     Purpose
 !     -------
-!     Interpolate K-coefficients to the given P,T and Qvar values.
+!     Interpolate K-coefficients to the given P,T.
 !
 !     Notes
@@ -58,8 +57,7 @@
       real*8 Tref(L_NTREF)
       real*8 pref(L_PINT)
-      real*8 wrefvar(L_REFVAR)
 
-      integer MT, MP, N, M, NP, NVAR
-      real*8  PW, TW, Qvar, wratio
+      integer MT, MP, N, M, NP
+      real*8  PW, TW
       real*8  PWL, LCOEF(4), T, U
 
@@ -143,24 +141,4 @@
       LCOEF(4) = (1.0-T)*U
 
-!  Get the indicies for abundance of the varying species. There are 10 sets of 
-!  k-coefficients with differing amounts of variable vs. constant gas.
-
-      IF(QVAR.le.WREFVAR(1)) then
-        NVAR   = 1
-        WRATIO = 0.0D0      ! put all the weight on the first point
-      ELSEIF(QVAR.ge.WREFVAR(L_REFVAR)) then
-        NVAR   = L_REFVAR-1 ! TB16 in order to not oversize NVAr when doing
-                                  !NVAR+1
-        WRATIO = 1.00D0     ! put all the weight on the last point
-      ELSE
-        DO N=2,L_REFVAR
-          IF(QVAR.GE.WREFVAR(N-1) .and. QVAR.lt.WREFVAR(N)) then
-            NVAR   = N-1
-            WRATIO = (QVAR - WREFVAR(N-1))/(WREFVAR(N) - WREFVAR(N-1))
-            GOTO 30
-          END IF
-        END DO
-      END IF
-
    30 CONTINUE
 
Index: trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1647)
+++ trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1648)
@@ -19,8 +19,7 @@
       real varian      ! Characteristic variance of log-normal distribution
       real r3n_q     ! used to compute r0 from number and mass mixing ratio
-      real rho_dust     ! Mars dust density (kg.m-3)
       real ref_r0        ! for computing reff=ref_r0*r0 (in log.n. distribution)
 !$OMP THREADPRIVATE(noms,mmol,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor, &
-	!$OMP varian,r3n_q,rho_dust,rho_ice,ref_r0)
+	!$OMP varian,r3n_q,ref_r0)
 
 ! tracer indexes: these are initialized in initracer and should be 0 if the
