Index: trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F	(revision 3275)
@@ -776,10 +776,10 @@
       DO iq=1,nqtot
           txt=trim(tname(iq))//"_surf"
-          if (txt.eq."h2o_vap") then
-            ! There is no surface tracer for h2o_vap;
+          if (txt.eq."h2o_gas") then
+            ! There is no surface tracer for h2o_gas;
             ! "h2o_ice" should be loaded instead
             txt="h2o_ice_surf"
             write(*,*) 'lect_start_archive: loading surface tracer',
-     &                     ' h2o_ice instead of h2o_vap'
+     &                     ' h2o_ice instead of h2o_gas'
           endif
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90	(revision 3275)
@@ -94,5 +94,5 @@
 
       real CLFtot
-      integer igen_ice,igen_vap ! to store the index of generic tracer
+      integer igen_ice,igen_gas ! to store the index of generic tracer
       logical dummy_bool ! dummy boolean just in case we need one
       ! integer i_rgcs_ice(aerogeneric)
Index: trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90	(revision 3275)
@@ -7,5 +7,5 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
           albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
-          tsurf,fract,dist_star,aerosol,muvar,                 &
+          zzlay,tsurf,fract,dist_star,aerosol,muvar,                 &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
@@ -30,10 +30,12 @@
       use aeropacity_mod, only: aeropacity
       use aeroptproperties_mod, only: aeroptproperties
-      use tracer_h, only: constants_epsi_generic
-      use comcstfi_mod, only: pi, mugaz, cpp
+      use tracer_h, only: constants_epsi_generic,igcm_ch4_gas,igcm_n2,mmol
+      use comcstfi_mod, only: pi, mugaz, cpp, r, g
       use callkeys_mod, only: varactive,diurnal,tracer,varfixed,satval, &
                               diagdtau,kastprof,strictboundcorrk,specOLR, &
                               tplanckmin,tplanckmax,global1d, &
-                              generic_condensation, aerohaze, haze_radproffix
+                              generic_condensation,aerohaze,haze_radproffix,&
+                              methane,carbox,cooling,nlte,strobel,&
+                              ch4fix,vmrch4_proffix,vmrch4fix
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -44,4 +46,5 @@
       use generic_tracer_index_mod, only: generic_tracer_index
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
+      use radcommon_h, only: wavev,wavei
       implicit none
 
@@ -80,4 +83,5 @@
       REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
       REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
+      REAL,INTENT(IN) :: zzlay(ngrid,nlayer)       ! Mid-layer altitude
       REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
       REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
@@ -157,8 +161,11 @@
 !$OMP THREADPRIVATE(tauaero)
       REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
+      real*8 nfluxtopv_nu(L_NSPECTV)
       REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
       REAL*8 nfluxtopi_nu(L_NSPECTI)                 ! Net band-resolved IR flux at TOA (W/m2).
       REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)         ! For 1D diagnostic.
       REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
+      real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) !
+      real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) !
       REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
       REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
@@ -211,9 +218,23 @@
       REAL*8 surface_stellar_flux   ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
 
+      !     NLTE factor for CH4
+      real eps_nlte_sw23(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
+      real eps_nlte_sw33(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
+      real eps_nlte_lw(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
+      integer Nfine,ifine
+      parameter(Nfine=701)
+      real,save :: levdat(Nfine),vmrdat(Nfine)
+      real :: vmrch4(ngrid,nlayer)              ! vmr ch4 from vmrch4_proffix
+
+      REAL dtlw_nu(nlayer,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands
+      REAL dtsw_nu(nlayer,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands
+
       ! local variable
+      REAL dpp  ! intermediate
+
       integer ok ! status (returned by NetCDF functions)
 
-      integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
-      logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
+      integer igcm_generic_gas, igcm_generic_ice! index of the vap and ice of generic_tracer
+      logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer
       real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
 !$OMP THREADPRIVATE(metallicity)
@@ -327,7 +348,5 @@
          ENDIF
 
-#ifndef MESOSCALE
          if (is_master) call system('rm -f surf_vals_long.out')
-#endif
 
          call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
@@ -434,5 +453,5 @@
 
 
-         if(varfixed .and. generic_condensation)then
+         if(varfixed .and. (generic_condensation .or. methane .or. carbox))then
             write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
             qvap_deep=-1. ! default value
@@ -504,4 +523,40 @@
       endif
 
+
+!-----------------------------------------------------------------------
+!     Prepare CH4 mixing ratio for radiative transfer
+      IF (methane) then
+        vmrch4(:,:)=0.
+
+        if (ch4fix) then
+           if (vmrch4_proffix) then
+            !! Interpolate on the model vertical grid
+             do ig=1,ngrid
+               CALL interp_line(levdat,vmrdat,Nfine, &
+                      zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
+             enddo
+           else
+            vmrch4(:,:)=vmrch4fix
+           endif
+        else
+            vmrch4(:,:)=pq(:,:,igcm_ch4_gas)*100.* &
+                        mmol(igcm_n2)/mmol(igcm_ch4_gas)
+        endif
+      ENDIF
+
+!     Prepare NON LTE correction in Pluto atmosphere
+      IF (nlte) then
+        CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,&
+                  eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
+      ENDIF
+!     Net atmospheric radiative cooling rate from C2H2 (K.s-1):
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    !   dtlw_hcn_c2h2=0.
+      if (cooling) then
+        !  CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,&
+                                !   pt,dtlw_hcn_c2h2)
+      endif
+
+
 !-----------------------------------------------------------------------
       do ig=1,ngrid ! Starting Big Loop over every GCM column
@@ -671,5 +726,5 @@
       !-----------------------------------------------------------------------
 
-      if (generic_condensation) then
+      if (generic_condensation .or. methane .or. carbox) then
 
          ! For now, only one GCS tracer can be both variable and radiatively active
@@ -679,11 +734,11 @@
          do iq=1,nq
 
-            call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-
-            if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
+            call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
+
+            if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
 
                if(varactive)then
 
-                  i_var=igcm_generic_vap
+                  i_var=igcm_generic_gas
                   do l=1,nlayer
                      qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
@@ -731,5 +786,5 @@
       !-----------------------------------------------------------------------
 
-      if (.not. generic_condensation) then
+      if (.not. (generic_condensation .or. methane .or. carbox)) then
          do k=1,L_LEVELS
             qvar(k) = 1.0D-7
@@ -741,12 +796,14 @@
          ! IMPORTANT: Now convert from kg/kg to mol/mol.
          do k=1,L_LEVELS
-            if (generic_condensation) then
+            if (generic_condensation .or. methane .or. carbox) then
                do iq=1,nq
-                  call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-                  if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
-
-                     epsi_generic=constants_epsi_generic(iq)
-
-                     qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic))
+                  call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
+                  if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
+                    if(.not. varactive .or. i_var.eq.iq)then
+
+                        epsi_generic=constants_epsi_generic(iq)
+
+                        qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic))
+                    endif
 
                   endif
@@ -762,5 +819,5 @@
       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
+             muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
       END DO
 
@@ -965,8 +1022,8 @@
                  tmid,pmid,taugsurf,qvar,muvarrad)
 
-            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
-                 acosz,stel_fract,                                 &
-                 nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,   &
-                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
+            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,   &
+                 acosz,stel_fract,nfluxtopv,fluxtopvdn,nfluxoutv_nu,&
+                 nfluxgndv_nu,nfluxtopv_nu, &
+                 fmnetv,fmnetv_nu,fluxupv,fluxdnv,fzerov,taugsurf)
 
          else ! During the night, fluxes = 0.
@@ -1010,5 +1067,5 @@
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
               wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu,         &
-              fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
+              fmneti,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
 
 !-----------------------------------------------------------------------
@@ -1050,17 +1107,35 @@
 
 !     Finally, the heating rates
-
          DO l=2,L_NLAYRAD
-            dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
-                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
-            dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
-                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
+            ! dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
+            !     *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
+            dpp = glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
+            do nw=1,L_NSPECTV
+                dtsw_nu(L_NLAYRAD+1-l,nw)= &
+                   (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp
+            end do
+
+            ! dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
+            !     *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
+            do nw=1,L_NSPECTI
+                dtlw_nu(L_NLAYRAD+1-l,nw)= &
+                   (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp
+            end do
          END DO
 
 !     These are values at top of atmosphere
-         dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
-             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
-         dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
-             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
+        !  dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
+        !      *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
+        !  dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
+        !      *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
+         dpp = g/(cpp*scalep*(plevrad(3)-plevrad(1)))
+         do nw=1,L_NSPECTV
+            dtsw_nu(L_NLAYRAD,nw)= &
+             (fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp
+         end do
+         do nw=1,L_NSPECTI
+             dtlw_nu(L_NLAYRAD,nw)= &
+             (fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp
+         end do
 
       !  Optical thickness diagnostics (added by JVO)
@@ -1084,9 +1159,67 @@
       endif
 
+! **********************************************************
+!     NON NLTE correction in Pluto atmosphere
+!     And conversion of LW spectral heating rates to total rates
+! **********************************************************
+
+      if (.not.nlte) then
+        eps_nlte_sw23(ig,:) =1. ! IF no NLTE
+        eps_nlte_sw33(ig,:) =1. ! IF no NLTE
+        eps_nlte_lw(ig,:) =1. ! IF no NLTE
+     endif
+
+     do l=1,nlayer
+
+        !LW
+        dtlw(ig,l) =0.
+!           dtlw_co(ig,l) =0.  ! only for diagnostic
+        do nw=1,L_NSPECTI
+          ! wewei : wavelength in micrometer
+          if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then
+            dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l)
+          else
+            !dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996)
+            dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996)
+!               dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic
+          end if
+          dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength
+        end do
+        ! adding c2h2 if cooling active
+        ! dtlw(ig,l)=dtlw(ig,l)+dtlw_hcn_c2h2(ig,l)
+
+        !SW
+        dtsw(ig,l) =0.
+
+        if (strobel) then
+
+         do nw=1,L_NSPECTV
+          if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then
+            dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
+          elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then
+            dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l)
+          else
+            dtsw_nu(l,nw)=dtsw_nu(l,nw)
+          end if
+          dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
+         end do
+
+        else ! total heating rates multiplied by nlte coef
+
+         do nw=1,L_NSPECTV
+            dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
+            dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
+         enddo
+
+        endif
+
+
+     end do
+! **********************************************************
+
 
 !-----------------------------------------------------------------------
       end do ! End of big loop over every GCM column.
 !-----------------------------------------------------------------------
-
 
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/condensation_generic_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/condensation_generic_mod.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/condensation_generic_mod.F90	(revision 3275)
@@ -18,5 +18,5 @@
 !     -------
 !     Calculates large-scale condensation of generic tracer "tname".
-!     By convention, tname ends with the suffix "_vap", as it represents the 
+!     By convention, tname ends with the suffix "_gas", as it represents the 
 !     gas phase of the generic tracer
 !     
@@ -39,5 +39,5 @@
         REAL, intent(in) :: pdt(ngrid,nlayer)     ! physical temperature tendency (K/s)
         REAL, intent(in) :: pdq(ngrid,nlayer,nq)  ! physical tracer tendency (K/s)
-        ! CHARACTER(*), intent(in) :: tname_vap     ! name of the tracer we consider. BY convention, it ends with _vap !!!
+        ! CHARACTER(*), intent(in) :: tname_gas     ! name of the tracer we consider. BY convention, it ends with _gas !!!
         REAL, intent(out) :: pdtlsc(ngrid,nlayer)  ! incrementation de la temperature (K)
         REAL, intent(out) :: pdqvaplsc(ngrid,nlayer,nq) ! incrementation de la vapeur du traceur
@@ -49,11 +49,11 @@
         REAL, SAVE :: qvap_deep   ! deep mixing ratio of vapor when simulating bottom less planets
         REAL, SAVE :: qvap_top   ! top mixing ratio of vapor when simulating bottom less planets
-        logical, save :: perfect_vap_profile
-!$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_vap_profile)
+        logical, save :: perfect_gas_profile
+!$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_gas_profile)
 
 !       Local variables
 
         ! to call only one time the ice/vap pair of a tracer
-        logical call_ice_vap_generic
+        logical call_ice_gas_generic
 
         INTEGER i, k , nn, iq
@@ -69,5 +69,5 @@
         real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
         real zqs_temp_1, zqs_temp_2, zqs_temp_3
-        integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
+        integer igcm_generic_gas, igcm_generic_ice! index of the vap and ice of generic_tracer
         ! CHARACTER(len=*) :: tname_ice
         ! evaporation calculations
@@ -96,8 +96,8 @@
                 write(*,*) " qvap_top = ",qvap_top 
                 
-                write(*,*) " perfect_vap_profile ? "
-                perfect_vap_profile=.false. ! default value
-                call getin_p("perfect_vap_profile",perfect_vap_profile)
-                write(*,*) " perfect_vap_profile = ",perfect_vap_profile
+                write(*,*) " perfect_gas_profile ? "
+                perfect_gas_profile=.false. ! default value
+                call getin_p("perfect_gas_profile",perfect_gas_profile)
+                write(*,*) " perfect_gas_profile = ",perfect_gas_profile
 
                 firstcall = .false.
@@ -115,9 +115,9 @@
         do iq=1,nq
 
-                call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-
-                if(call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
+                call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
+
+                if(call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
                         m=constants_mass(iq)
-                        delta_vapH=constants_delta_vapH(iq)
+                        delta_gasH=constants_delta_gasH(iq)
                         Tref=constants_Tref(iq)
                         Pref=constants_Pref(iq)
@@ -139,5 +139,5 @@
                                         !	    zt(i)=15.   ! check too low temperatures
                                         endif
-                                        zx_q(i) = pq(i,k,igcm_generic_vap)+pdq(i,k,igcm_generic_vap)*ptimestep
+                                        zx_q(i) = pq(i,k,igcm_generic_gas)+pdq(i,k,igcm_generic_gas)*ptimestep
 
                                         call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
@@ -189,12 +189,12 @@
 
                                 !Tendances de t et q
-                                pdqvaplsc(1:ngrid,k,igcm_generic_vap)  = - zcond(1:ngrid) 
-                                pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_vap)
+                                pdqvaplsc(1:ngrid,k,igcm_generic_gas)  = - zcond(1:ngrid) 
+                                pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_gas)
                                 pdtlsc(1:ngrid,k)  = pdtlsc(1:ngrid,k) + pdqliqlsc(1:ngrid,k,igcm_generic_ice)*Lcp
 
                         Enddo ! k= nlayer, 1, -1
 
-                        if ((perfect_vap_profile) .and. (ngrid.eq.1)) then
-                                ! perfect_vap_profile is a mode that should a priori only be used in 1D:
+                        if ((perfect_gas_profile) .and. (ngrid.eq.1)) then
+                                ! perfect_gas_profile is a mode that should a priori only be used in 1D:
                                 ! as it is written below, it aims to force the vap profile:
                                 ! - below condensation, it is forced to qvap_deep
@@ -202,5 +202,5 @@
                                 ! - above the cold trap, it is forced to the value allowed by the cold trap
 
-                                ! perfect_vap_profile can be customed as you want
+                                ! perfect_gas_profile can be customed as you want
                                 
                                 tau = 10. * ptimestep ! tau must not be lower than the physical step time. 
@@ -223,5 +223,5 @@
                                 if (k_cold_trap .lt. nlayer) then
                                         do k=k_cold_trap+1,nlayer
-                                                pdqvaplsc(1,k,igcm_generic_vap) = (pq(1,k_cold_trap,igcm_generic_vap) - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
+                                                pdqvaplsc(1,k,igcm_generic_gas) = (pq(1,k_cold_trap,igcm_generic_gas) - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
                                         enddo
                                 endif
@@ -231,8 +231,8 @@
                                         call Psat_generic(zt(1),pplay(1,k),metallicity,psat_tmp,zqs_temp)
                                         if (zqs_temp .gt. qvap_deep) then
-                                                pdqvaplsc(1,k,igcm_generic_vap)  = (qvap_deep - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
+                                                pdqvaplsc(1,k,igcm_generic_gas)  = (qvap_deep - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
                                         endif
                                         if (zqs_temp .lt. qvap_deep) then
-                                                pdqvaplsc(1,k,igcm_generic_vap)  = (0.99*zqs_temp - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)
+                                                pdqvaplsc(1,k,igcm_generic_gas)  = (0.99*zqs_temp - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
                                         endif
                                 enddo
@@ -246,5 +246,5 @@
                                 ! brings lower generic vapor ratio to a fixed value.
                                 ! tau is in seconds and must not be lower than the physical step time.
-                                pdqvaplsc(1:ngrid,1,igcm_generic_vap) = (qvap_deep - pq(1:ngrid,1,igcm_generic_vap))/tau - pdq(1:ngrid,1,igcm_generic_vap)
+                                pdqvaplsc(1:ngrid,1,igcm_generic_gas) = (qvap_deep - pq(1:ngrid,1,igcm_generic_gas))/tau - pdq(1:ngrid,1,igcm_generic_gas)
                         endif
                         if (qvap_top >= 0.) then
@@ -252,7 +252,7 @@
                                 ! brings lower generic vapor ratio to a fixed value.
                                 ! tau is in seconds and must not be lower than the physical step time.
-                                pdqvaplsc(1:ngrid,nlayer,igcm_generic_vap) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_vap))/tau - pdq(1:ngrid,nlayer,igcm_generic_vap)
+                                pdqvaplsc(1:ngrid,nlayer,igcm_generic_gas) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_gas))/tau - pdq(1:ngrid,nlayer,igcm_generic_gas)
                         endif
-                endif !if(call_ice_vap_generic)
+                endif !if(call_ice_gas_generic)
         enddo ! iq=1,nq
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/convadj.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/convadj.F	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/convadj.F	(revision 3275)
@@ -304,10 +304,10 @@
          end do
          ! do l = 1, nlay
-         !    print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_vap)
+         !    print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_gas)
          ! end do
          ! do l = 1, nlay
-         !    print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_vap)
+         !    print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_gas)
          ! end do
-         !    print*,'zqm(vap)         = ',zqm(igcm_h2o_vap)
+         !    print*,'zqm(vap)         = ',zqm(igcm_h2o_gas)
             print*,'jadrs=',jadrs
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/generic_cloud_common_h.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/generic_cloud_common_h.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/generic_cloud_common_h.F90	(revision 3275)
@@ -16,5 +16,5 @@
     implicit none 
     real, save :: m  ! molecular mass of the specie (g/mol)
-    real, save :: delta_vapH ! Enthalpy of vaporization (J/mol)
+    real, save :: delta_gasH ! Enthalpy of vaporization (J/mol)
     real,save :: Tref ! Ref temperature for Clausis-Clapeyron (K)
     real,save :: Pref ! Reference pressure for Clausius Clapeyron (Pa)
@@ -46,62 +46,62 @@
         if (trim(specname) .eq. "Mg") then
             print*,"Loading data for Mg"
-            delta_vapH = 521.7E3 !J/mol
+            delta_gasH = 521.7E3 !J/mol
             Tref = 2303.0 !K
             Pref =   1.0E5 !in Pa
             m = 140.6931
-            RLVTT_generic = delta_vapH/(m/1000.)
+            RLVTT_generic = delta_gasH/(m/1000.)
             metallicity_coeff=1.0*log(10.0)
             
         else if (trim(specname) .eq. "Na") then
             print*,"Loading data for Na"
-            delta_vapH = 265.9E3
+            delta_gasH = 265.9E3
             Tref = 1624.0
             Pref =  1.0E5
             m = 78.04454 !Na2S
-            RLVTT_generic = delta_vapH/(m/1000.)
+            RLVTT_generic = delta_gasH/(m/1000.)
             metallicity_coeff=0.5*log(10.0)
             
         else if (trim(specname) .eq. "Fe") then 
             print*,"Loading data for Fe"
-            delta_vapH = 401.9E3 
+            delta_gasH = 401.9E3 
             Tref = 2903.0
             Pref  = 1.0E5
             m = 55.8450
-            RLVTT_generic = delta_vapH/(m/1000.)
+            RLVTT_generic = delta_gasH/(m/1000.)
             metallicity_coeff=0.0*log(10.0)
             
         else if (trim(specname) .eq. "Cr") then
             print*,"Loading data for Cr"
-            delta_vapH = 394.2E3
+            delta_gasH = 394.2E3
             Tref = 2749.0
             Pref  = 1.0E5
             m = 51.9961
-            RLVTT_generic = delta_vapH/(m/1000.)
+            RLVTT_generic = delta_gasH/(m/1000.)
             metallicity_coeff=0.0*log(10.0)
             
         else if (trim(specname) .eq. "KCl") then 
             print*,"Loading data for KCl"
-            delta_vapH = 217.9E3
+            delta_gasH = 217.9E3
             Tref = 1495.0
             Pref = 1.0E5
             metallicity_coeff=0.0*log(10.0)
             m = 74.5498
-            RLVTT_generic = delta_vapH/(m/1000.)
+            RLVTT_generic = delta_gasH/(m/1000.)
         else if (trim(specname) .eq. "Mn") then 
             print*,"Loading data for Mn"
-            delta_vapH = 455.8E3
+            delta_gasH = 455.8E3
             Tref = 2064.0
             Pref = 1.0E5
             metallicity_coeff=1.0*log(10.0)
             m = 87.003049
-            RLVTT_generic = delta_vapH/(m/1000.)
+            RLVTT_generic = delta_gasH/(m/1000.)
         else if (trim(specname) .eq. "Zn") then 
             print*,"Loading data for Zn"
-            delta_vapH = 303.9E3
+            delta_gasH = 303.9E3
             Tref = 1238.0
             Pref = 1.0E5
             metallicity_coeff=1.0*log(10.0)
             m = 97.445
-            RLVTT_generic = delta_vapH/(m/1000.)
+            RLVTT_generic = delta_gasH/(m/1000.)
         else
             print*,"Unknow species (not in Mg, Fe, Na, KCl, Cr, Mn or Zn)"
@@ -138,7 +138,7 @@
             if (index(table_traceurs_line,specname) /= 0) then
                 
-                if (index(table_traceurs_line,'delta_vapH='   ) /= 0) then
-                    read(table_traceurs_line(index(table_traceurs_line,'delta_vapH=')+len('delta_vapH='):),*) delta_vapH
-                    print*, 'delta_vapH ', delta_vapH
+                if (index(table_traceurs_line,'delta_gasH='   ) /= 0) then
+                    read(table_traceurs_line(index(table_traceurs_line,'delta_gasH=')+len('delta_gasH='):),*) delta_gasH
+                    print*, 'delta_gasH ', delta_gasH
                 end if
                 if (index(table_traceurs_line,'Tref='   ) /= 0) then
@@ -178,5 +178,5 @@
         close(117)
 
-        RLVTT_generic=delta_vapH/(m/1000.)
+        RLVTT_generic=delta_gasH/(m/1000.)
 
         write(*,*) 'RLVTT_generic', RLVTT_generic
@@ -190,6 +190,6 @@
     !============================================================================
     !   Clausius-Clapeyron relation : 
-    !   d(Ln(psat))/dT = delta_vapH/(RT^2)
-    !   ->Psat = Pref * exp(-delta_vapH/R*(1/T-1/Tref))
+    !   d(Ln(psat))/dT = delta_gasH/(RT^2)
+    !   ->Psat = Pref * exp(-delta_gasH/R*(1/T-1/Tref))
     !
     !   Authors 
@@ -203,5 +203,5 @@
         real, intent(out) :: psat,qsat !saturation vapor pressure (Pa) and mass mixing ratio at saturation (kg/kg) of the layer 
 
-        psat = pref*exp(-delta_vapH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)
+        psat = pref*exp(-delta_gasH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)
 
         if (psat .gt. p) then 
@@ -217,5 +217,5 @@
     !===============================================================================
     !   Compute dqsat = L/cp* d(qsat)/dT and d(Ln(psat)) = L/cp*d(Ln(psat))/dT
-    !   we have d(ln(psat))/dT =  delta_vapH/R*(1/T^2) for clausius-clapeyron
+    !   we have d(ln(psat))/dT =  delta_gasH/R*(1/T^2) for clausius-clapeyron
     !   r*mugaz/1000. is the perfect gaz constant in the computation of "dummy"
     !   Authors 
@@ -235,5 +235,5 @@
         real :: dummy ! used to store d(ln(psat))/dT
 
-        dummy = delta_vapH/((r*mugaz/1000.)*(T**2)) 
+        dummy = delta_gasH/((r*mugaz/1000.)*(T**2)) 
 
         if (psat .gt. p) then 
@@ -270,5 +270,5 @@
 
 
-        Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_vapH*Log(p/Pref))
+        Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_gasH*Log(p/Pref))
 
     end subroutine Tsat_generic
Index: trunk/LMDZ.PLUTO/libf/phypluto/generic_tracer_index_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/generic_tracer_index_mod.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/generic_tracer_index_mod.F90	(revision 3275)
@@ -6,6 +6,6 @@
 !     for both the vap & ice parts of the tracer
 !
-!     call_ice_vap_generic is a boolean which tells if the tracer is condensable
-!     and allows, if it is condensable, to call the vap/ice pair 
+!     call_ice_gas_generic is a boolean which tells if the tracer is condensable
+!     and allows, if it is condensable, to call the vap/ice pair
 !
 !     Noé Clément & Lucas Teinturier (2022)
@@ -17,15 +17,15 @@
     contains
 
-    subroutine generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
+    subroutine generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
 
         USE tracer_h, only: noms, is_condensable
-        integer nq,iq, igcm_generic_vap, igcm_generic_ice, i_shortname, ii
-        logical call_ice_vap_generic
+        integer nq,iq, igcm_generic_gas, igcm_generic_ice, i_shortname, ii
+        logical call_ice_gas_generic
         character(len=10) :: shortname
-        if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0) &
+        if((is_condensable(iq)==1) .and. (index(noms(iq),"gas") .ne. 0) &
                                                     .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
-            ! Let's get the index of our tracers (we look for igcm _generic_vap and igcm_generic_ice)
+            ! Let's get the index of our tracers (we look for igcm _generic_gas and igcm_generic_ice)
 
-            igcm_generic_vap = iq
+            igcm_generic_gas = iq
             igcm_generic_ice = -1
 
@@ -46,5 +46,5 @@
             ! ! endif
 
-            ! ! call_ice_vap_generic = .true.
+            ! ! call_ice_gas_generic = .true.
             i_shortname = index(noms(iq),'_')-1
             shortname=noms(iq)(1:i_shortname)
@@ -59,7 +59,7 @@
             endif
             ! if we didn't find it before, we look after the vap tracer in pq
-            if (igcm_generic_ice .eq. -1 .and. (iq<nq)) then 
+            if (igcm_generic_ice .eq. -1 .and. (iq<nq)) then
                 do ii=iq+1,nq
-                    if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then 
+                    if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then
                         igcm_generic_ice=ii
                         exit
@@ -67,8 +67,8 @@
                 enddo
             endif
-            call_ice_vap_generic = .true.
+            call_ice_gas_generic = .true.
         else
 
-            call_ice_vap_generic = .false.
+            call_ice_gas_generic = .false.
 
         end if ! if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))
Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90	(revision 3275)
@@ -883,4 +883,23 @@
       " albmin_ch4 = ",albmin_ch4
 
+
+     if (is_master)write(*,*)trim(rname)//&
+     "fixed gaseous CH4 mixing ratio for rad transfer ?"
+     ch4fix=.false. ! default value
+     call getin_p("ch4fix",ch4fix)
+     if (is_master)write(*,*)trim(rname)//&
+       " ch4fix = ",ch4fix
+     if (is_master)write(*,*)trim(rname)//&
+       "fixed gaseous CH4 mixing ratio for rad transfer ?"
+     vmrch4fix=0.5 ! default value
+     call getin_p("vmrch4fix",vmrch4fix)
+     if (is_master)write(*,*)trim(rname)//&
+        " vmrch4fix = ",vmrch4fix
+     vmrch4_proffix=.false. ! default value
+     call getin_p("vmrch4_proffix",vmrch4_proffix)
+     if (is_master)write(*,*)trim(rname)//&
+        " vmrch4_proffix = ",vmrch4_proffix
+
+
 !*********** CO *********************************
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90	(revision 3275)
@@ -127,5 +127,5 @@
        IF (.NOT. allocated(is_rgcs)) allocate(is_rgcs(nq)) !LT
        IF (.NOT. allocated(constants_mass)) allocate(constants_mass(nq))
-       IF (.NOT. allocated(constants_delta_vapH)) allocate(constants_delta_vapH(nq))
+       IF (.NOT. allocated(constants_delta_gasH)) allocate(constants_delta_gasH(nq))
        IF (.NOT. allocated(constants_Tref)) allocate(constants_Tref(nq))
        IF (.NOT. allocated(constants_Pref)) allocate(constants_Pref(nq))
@@ -157,5 +157,5 @@
        is_rgcs(:) = 0
        constants_mass(:)=0
-       constants_delta_vapH(:)=0
+       constants_delta_gasH(:)=0
        constants_Tref(:)=0
        constants_Pref(:)=0
@@ -229,5 +229,5 @@
           write(*,*) 'Tracer ',count,' = n2'
         endif
-        if (noms(iq).eq."ch4_vap") then
+        if (noms(iq).eq."ch4_gas") then
           igcm_ch4_gas=iq
           mmol(igcm_ch4_gas)=16.
@@ -243,5 +243,5 @@
           write(*,*) 'Tracer ',count,' = ch4 ice'
         endif
-        if (noms(iq).eq."co_vap") then
+        if (noms(iq).eq."co_gas") then
           igcm_co_gas=iq
           mmol(igcm_co_gas)=28.
@@ -376,5 +376,5 @@
                 call specie_parameters_table(noms(iq)(1:len(trim(noms(iq)))-4))
                 constants_mass(iq)=m
-                constants_delta_vapH(iq)=delta_vapH
+                constants_delta_gasH(iq)=delta_gasH
                 constants_Tref(iq)=Tref
                 constants_Pref(iq)=Pref
Index: trunk/LMDZ.PLUTO/libf/phypluto/mass_redistribution_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/mass_redistribution_mod.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/mass_redistribution_mod.F90	(revision 3275)
@@ -7,9 +7,9 @@
 		       pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr,  &
 		       pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr)
-                                                   
+
        use radcommon_h, only: glat
        USE planete_mod, only: bp
        use comcstfi_mod, only: g
-       
+
        IMPLICIT NONE
 !=======================================================================
@@ -22,11 +22,11 @@
 !
 !   input:
-!   ----- 
+!   -----
 !    ngrid                 nombre de points de verticales
 !                          (toutes les boucles de la physique sont au
 !                          moins vectorisees sur ngrid)
 !    nlayer                nombre de couches
-!    pplay(ngrid,nlayer)   Pressure levels 
-!    pplev(ngrid,nlayer+1) Pressure levels 
+!    pplay(ngrid,nlayer)   Pressure levels
+!    pplev(ngrid,nlayer+1) Pressure levels
 !    nq                    Number of tracers
 !
@@ -35,5 +35,5 @@
 !    pu,pv (ngrid,nlayer)  wind velocity (m/s)
 !
-!                    
+!
 !    pdX                   physical tendency of X before mass redistribution
 !
@@ -45,5 +45,5 @@
 !    pdXmr(ngrid)           physical tendency of X after mass redistribution
 !
-!    
+!
 !
 !=======================================================================
@@ -55,5 +55,5 @@
 !    Arguments :
 !    ---------
-      INTEGER,INTENT(IN) :: ngrid, nlayer, nq   
+      INTEGER,INTENT(IN) :: ngrid, nlayer, nq
       REAL,INTENT(IN) :: ptimestep
       REAL,INTENT(IN) :: pcapcal(ngrid)
@@ -80,5 +80,5 @@
       REAL zzu(nlayer),zzv(nlayer)
       REAL zzq(nlayer,nq),zzt(nlayer)
-!    Dummy variables      
+!    Dummy variables
       INTEGER n,l,ig,iq
       REAL zdtsig(ngrid,nlayer)
@@ -88,5 +88,5 @@
       REAL zq1(nlayer)
       REAL ztsrf(ngrid)
-      REAL ztm(nlayer+1) 
+      REAL ztm(nlayer+1)
       REAL zum(nlayer+1) , zvm(nlayer+1)
       REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
@@ -103,11 +103,11 @@
 !
       IF (firstcall) THEN
-         firstcall=.false.	 
+         firstcall=.false.
       ENDIF
 !
 !======================================================================
-!    Calcul of h2o condensation  
+!    Calcul of h2o condensation
 !    ============================================================
-!  
+!
 !    Used variable :
 !       pdmassmr      : air Mass added to the atmosphere in each layer per unit time (kg/m2/s)
@@ -135,5 +135,8 @@
 !    Changing pressure at the surface:
 !    """"""""""""""""""""""""""""""""""""
-         
+
+      zmassboil(1:ngrid) = 0 ! AF: TODO: update this for pluto
+      pdtsrfmr(1:ngrid) = 0
+
       pdpsrfmr(1:ngrid) = (zdmass_sum(1:ngrid,1)+zmassboil(1:ngrid))*g
 
@@ -151,5 +154,5 @@
 
 ! ***************************************************************
-!  Correction to account for redistribution between sigma or hybrid 
+!  Correction to account for redistribution between sigma or hybrid
 !  layers when changing surface pressure
 !  zzx quantities have dimension (nlayer) to speed up calculation
@@ -174,10 +177,10 @@
            ! Ehouarn: but for now leave things as before
             zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1))
-! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
+! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
             if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
          END DO
 
 ! Mass of each layer
-! ------------------ 
+! ------------------
          zzmass(1:nlayer)=zmass(ig,1:nlayer)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1))
 
@@ -186,19 +189,19 @@
 !  """"""""""""""""""""""""""""""""
 
-!        averaging operator for TRANSPORT  
+!        averaging operator for TRANSPORT
 !        """"""""""""""""""""""""""""""""
 !        Value transfert at the surface interface when condensation
 !        sublimation:
          ztm(1) = ztsrf(ig)
-         zum(1) = 0.  
-         zvm(1) = 0.  
+         zum(1) = 0.
+         zvm(1) = 0.
          zqm(1,1:nq)=0. ! most tracer do not condense !
-         !! if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
-	 
+         !! if (water) zqm(1,igcm_h2o_gas)=1. ! flux is 100% h2o at surface
+
 !        Van Leer scheme:
          w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep
-         call vl1d(nlayer,zzt,2.,zzmass,w,ztm) 
-         call vl1d(nlayer,zzu,2.,zzmass,w,zum) 
-         call vl1d(nlayer,zzv,2.,zzmass,w,zvm) 
+         call vl1d(nlayer,zzt,2.,zzmass,w,ztm)
+         call vl1d(nlayer,zzu,2.,zzmass,w,zum)
+         call vl1d(nlayer,zzv,2.,zzmass,w,zvm)
          do iq=1,nq
            zq1(1:nlayer)=zzq(1:nlayer,iq)
@@ -214,5 +217,5 @@
 
 !        Surface condensation affects low winds
-         if (zmflux(1).lt.0) then 
+         if (zmflux(1).lt.0) then
             zum(1)= zzu(1) *  (w(1)/zzmass(1))
             zvm(1)= zzv(1) *  (w(1)/zzmass(1))
@@ -222,14 +225,14 @@
             end if
          end if
-                    
-         ztm(nlayer+1)= zzt(nlayer) ! should not be used, but... 
+
+         ztm(nlayer+1)= zzt(nlayer) ! should not be used, but...
          zum(nlayer+1)= zzu(nlayer)  ! should not be used, but...
          zvm(nlayer+1)= zzv(nlayer)  ! should not be used, but...
          zqm(nlayer+1,1:nq)= zzq(nlayer,1:nq)
- 
-!        Tendencies on T, U, V, Q 
+
+!        Tendencies on T, U, V, Q
 !        """"""""""""""""""""""""
          DO l=1,nlayer
- 
+
 !           Tendencies on T
             pdtmr(ig,l) = (1/zzmass(l)) *   &
@@ -253,5 +256,5 @@
          enddo
 
-      END DO  ! loop on ig 
+      END DO  ! loop on ig
 
 CONTAINS
@@ -260,5 +263,5 @@
       SUBROUTINE vl1d(llm,q,pente_max,zzmass,w,qm)
 !
-!    
+!
 !     Operateur de moyenne inter-couche pour calcul de transport type
 !     Van-Leer " pseudo amont " dans la verticale
@@ -270,5 +273,5 @@
 !
 !   --------------------------------------------------------------------
-      
+
       IMPLICIT NONE
 
@@ -280,5 +283,5 @@
       REAL w(llm+1)
 !
-!      Local 
+!      Local
 !   ---------
 !
@@ -287,9 +290,9 @@
       real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
       real sigw, Mtot, MQtot
-      integer m 
-!     integer ismax,ismin 
-
-
-!    On oriente tout dans le sens de la pression 
+      integer m
+!     integer ismax,ismin
+
+
+!    On oriente tout dans le sens de la pression
 !     W > 0 WHEN DOWN !!!!!!!!!!!!!
 
@@ -347,6 +350,6 @@
                stop
              end if
-          else      ! if(w(l+1).lt.0) 
-             m = l-1 
+          else      ! if(w(l+1).lt.0)
+             m = l-1
              Mtot = zzmass(m+1)
              MQtot = zzmass(m+1)*q(m+1)
@@ -365,5 +368,5 @@
              else
                 qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
-             end if   
+             end if
           end if
        enddo
@@ -373,5 +376,5 @@
 !         if(w(1).gt.0.) then
 !            qm(1)=q(1)
-!         else 
+!         else
 !           qm(1)=0.
 !         end if
Index: trunk/LMDZ.PLUTO/libf/phypluto/nlte_ch4.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/nlte_ch4.F90	(revision 3275)
+++ trunk/LMDZ.PLUTO/libf/phypluto/nlte_ch4.F90	(revision 3275)
@@ -0,0 +1,304 @@
+      subroutine nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,eps23,eps33,eps76)
+
+      use comgeomfi_h
+      use callkeys_mod, only: strobel
+      use tracer_h, only: igcm_n2,mmol
+      use comcstfi_mod, only: pi, g, cpp,r
+      implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Calculation of epsilon - non LTE efficiency factor
+!     For CH4 heating and cooling rate calculation
+!     Cf Strobel et al 1996
+!
+!     Inputs
+!     ------
+!     ngrid                 Number of vertical columns
+!     nlayer                Number of layers
+!     pplay(ngrid,nlayer)   Pressure layers
+!     pt
+!
+!     Outputs
+!     -------
+!
+!     eps_sw       ! efficiency factor for heating rate
+!     eps_lw       ! efficiency factor for cooling rate
+!
+!     Authors
+!     -------
+!     Tanguy Bertrand (2015,2022)
+!
+!==================================================================
+
+!-----------------------------------------------------------------------
+
+!     Arguments
+
+      INTEGER ngrid, nlayer, nq
+      REAL pplay(ngrid,nlayer)
+      REAL pplev(ngrid,nlayer)
+      REAL pt(ngrid,nlayer)
+      REAL vmrch4(ngrid,nlayer)
+
+!-----------------------------------------------------------------------
+!     Local variables
+
+      INTEGER l,ig
+
+      REAL cbol         ! Boltzman constant
+      REAL avog         ! avogadro constant
+      REAL beta23
+      REAL beta33
+      REAL nd(ngrid,nlayer)      ! density number  cm-3
+      REAL na(ngrid,nlayer)      ! column density number  cm-2
+      REAL rho(ngrid,nlayer)     ! density of atmosphere
+      REAL kn3,k3,k2n4,k3n4,k4
+      REAL An3,A3,A2n4,A3n4,A4
+      REAL dphin3,dphi3,dphi2n4,dphi3n4,dphi4
+      REAL phin3(ngrid,nlayer)
+      REAL phi3(ngrid,nlayer)
+      REAL phi2n4(ngrid,nlayer)
+      REAL phi3n4(ngrid,nlayer)
+      REAL phi4(ngrid,nlayer)
+      REAL phi4T(ngrid,nlayer)
+
+!     DATA from Strobel et al 1996
+      REAL, DIMENSION(22) :: gam2       ! band escape fonction 2.3
+      REAL, DIMENSION(22) :: Naval2     ! column density number 2.3
+      REAL, DIMENSION(22) :: gam3       ! band escape fonction 3.3
+      REAL, DIMENSION(22) :: Naval3    ! column density number 3.3
+      REAL gammafi23(nlayer),gammafi33(nlayer)       ! Band espace fonction
+      REAL gammafi23mod(nlayer),gammafi33mod(nlayer)       ! Band espace fonction
+      REAL vmr
+      REAL k4zT  ! function to calculate k4 depending on T and vmr
+      REAL log_square ! log function
+      REAL log_affine ! log function
+      REAL ratio_affine ! linear approx function to calculate ratio for gamma calculation
+      REAL ratio_square !
+      REAL ratio_sq,ratio_af
+
+!     Output
+      REAL,intent(out) :: eps23(ngrid,nlayer)
+      REAL,intent(out) :: eps33(ngrid,nlayer)
+      REAL,intent(out) :: eps76(ngrid,nlayer)
+
+!     Option simplified coefficient
+      REAL klw, ksw,  eps_sw1Pa
+
+!-----------------------------------------------------------------------
+      eps23(:,:)=1.
+      eps33(:,:)=1.
+      eps76(:,:)=1.
+
+      IF (strobel) then
+
+!       Constantes
+        cbol=1.3806488e-23 !Boltzman constant
+        avog=6.022141e23
+
+        ! band 3
+        k3=2.8e-11    !cm-3
+        A3=25.2
+        dphi3=k3/A3
+        ! band n3
+        kn3=2.8e-11    !cm-3
+        An3=27.3
+        dphin3=kn3/An3
+        ! band 3n4
+        k3n4=1.16e-11  !cm-3
+        A3n4=2.12
+        dphi3n4=k3n4/(3.*A3n4)
+        ! band 4
+        !! Here we use function k4zT instead of constant k4
+        k4=4.5e-16        !cm-3
+        A4=2.12
+        dphi4=k4/A4
+        ! band 2n4
+        k2n4=1.16e-11
+        A2n4=2.12
+        dphi2n4=k2n4/(2.*A2n4)
+
+!       2.3 band :
+        beta23=0.09
+!       3.3 band :
+        beta33=0.13
+
+!       Escape Function Gamma (Data Strobel et al) :
+        ! 100 K 2.3
+        Naval2=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, &
+                7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, &
+                9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, &
+                1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/)
+        gam2=(/1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,9.90000000e-01, &
+              9.80000000e-01,9.50000000e-0,9.00000000e-01,8.00000000e-01,7.00000000e-01,5.00000000e-01, &
+              3.50000000e-01,2.20000000e-01,1.20000000e-01,6.50000000e-02,3.20000000e-02,1.40000000e-02, &
+              6.20000000e-03,2.50000000e-03,1.05000000e-03,4.83227951e-06/)
+
+        ! 100 K 3.3
+        Naval3=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, &
+                7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, &
+                9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, &
+                1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/)
+        gam3=(/1.00000000e+00,9.80000000e-01,9.00000000e-01,8.50000000e-01,7.00000000e-01,5.00000000e-01, &
+              3.00000000e-01,1.50000000e-01,1.00000000e-01,6.00000000e-02,4.00000000e-02,2.30000000e-02, &
+              1.45000000e-02,9.00000000e-03,5.00000000e-03,3.00000000e-03,1.50000000e-03,8.00000000e-04, &
+              3.80000000e-04,1.80000000e-04,8.30000000e-05,6.81737544e-07/)
+
+        DO ig=1,ngrid
+
+          phi3n4(ig,:)=1.
+          phi2n4(ig,:)=1.
+          phin3(ig,:)=1.
+          phi4(ig,:)=1.
+          nd(ig,:)=1.
+          na(ig,:)=1.
+
+          !write(*,*) 'TB22 pplev=',pplev(1,:)
+          !write(*,*) 'TB22 tlev=',pt(1,:)
+
+          DO l=1,nlayer
+           !!calculation of density
+           rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))
+
+           !!calculation of column number density for ch4 (cm-2)
+           na(ig,l)=vmrch4(ig,l)/100./(mmol(igcm_n2)*1.e-3)* &
+                             pplev(ig,l)/g*avog*1.e-4
+
+           !!calculation of atmospheric number density (cm-3) for each layer
+           nd(ig,l)=rho(ig,l)/(mmol(igcm_n2)*1.e-3)*avog*1.e-6
+
+           !!calculation of the phi :
+           phin3(ig,l)=dphin3*nd(ig,l)
+           phi3(ig,l)=dphi3*nd(ig,l)
+           phi2n4(ig,l)=dphi2n4*nd(ig,l)
+           phi3n4(ig,l)=dphi3n4*nd(ig,l)
+           phi4(ig,l)=dphi4*nd(ig,l)
+           !! Getting phi4 depending on T and vmr
+           vmr=vmrch4(ig,l)/100.
+           phi4T(ig,l)=k4zT(pt(ig,l),vmr)/A4*nd(ig,l)
+          ENDDO
+          !write(*,*) 'TB22 nd=',nd(1,:)
+          !write(*,*) 'TB22 vmr=',vmrch4
+          !write(*,*) 'TB22 phi4T=',phi4T
+
+          !! interpolation of na values in cm-2 from Strobel values for
+          !the 100K model:
+          CALL interp_line(Naval2,gam2,22,na(ig,:),gammafi23,nlayer)
+          CALL interp_line(Naval3,gam3,22,na(ig,:),gammafi33,nlayer)
+
+          !write(*,*) 'TB22rad na = ',na(1,:)
+          !write(*,*) 'TB22rad gammafi = ',gammafi23
+
+          DO l=1,nlayer
+
+           !! Getting ratio depending on temperature and n (cm-2)
+           ratio_sq=ratio_square(pt(ig,l),na(ig,l))
+           ratio_af=ratio_affine(pt(ig,l),na(ig,l))
+           !! Modified gamma depending on T
+           gammafi23mod(l)=min(ratio_sq*gammafi23(l),1.)
+           gammafi33mod(l)=min(ratio_af*gammafi33(l),1.)
+
+           !! Heating : cf Strobel 1996 and Zalucha 2011
+           ! eps 2.3
+           eps23(ig,l)=phin3(ig,l)/(phin3(ig,l)+gammafi23mod(l))*     &
+                      (beta23+phi3n4(ig,l)/(phi3n4(ig,l)+1)*  &
+                      phi4T(ig,l)/(phi4T(ig,l)+gammafi23mod(l))*   &
+                      (3.*phi2n4(ig,l)+1.)/(1.+phi2n4(ig,l))*(1-beta23)/3.)
+
+           ! eps 3.3
+           eps33(ig,l)=phi3(ig,l)/(phi3(ig,l)+gammafi33mod(l))*     &
+                      (beta33+phi2n4(ig,l)/(phi2n4(ig,l)+1)*  &
+                      phi4T(ig,l)/(phi4T(ig,l)+gammafi33mod(l))*   &
+                      (1-beta33))
+
+           ! Cooling from Zalucha et al 2013 eq10
+           eps76(ig,l)=1./(1.+gammafi23mod(l)/(2.*phi4T(ig,l)))
+
+          ENDDO
+          !write(*,*) 'TB22rad eps23 = ',eps23(1,:)
+          !write(*,*) 'TB22rad eps33 = ',eps33(1,:)
+          !write(*,*) 'TB22rad gamma33 = ',gammafi33mod
+          !write(*,*) 'TB22rad eps76 = ',eps76(1,:)
+
+        ENDDO   ! ngrid
+
+      ELSE     !    if (strobel)
+         !if strobel = false
+
+         klw = 3.6E-4  ! fit strobel
+         !eps_sw1Pa = 0.6 ! fit strobel
+         eps_sw1Pa = 0.90 ! to fit NH temperature en 1D [ch4]=0.5
+         !eps_sw1Pa = 0.98 ! to fit NH temperature en 3D = fit 1D [ch4]=0.25
+
+         ksw = 8.5E-5*(0.9/(eps_sw1Pa-0.1)-1.)
+         do l=1, nlayer
+           do ig=1, ngrid
+              rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))
+              eps23(ig,l) =  0.1+ 0.9/(1.+ksw/rho(ig,l))
+              eps76(ig,l) = 1./(1.+klw/rho(ig,l))
+           end do
+         end do
+
+      END IF
+
+      end subroutine nlte_ch4
+
+      real function k4zT(temp,vmr)
+      implicit none
+      real, intent(in) ::  temp,vmr
+      k4zT=1.9e-15*exp((temp-240.)/70.)*(1+10.*vmr)
+      return
+      end function k4zT
+
+      !! Log log affine and square function for ratio adjustment and extrapolation
+
+      real function log_square(x,y1,y2,x1,x2,coef)
+      implicit none
+      real, intent(in) ::  x,y1,y2,x1,x2,coef
+      real aa,bb
+      aa=(log(y2)-log(y1))/((log(coef*x2))**4-(log(coef*x1)**4))
+      bb=log(y2)-aa*(log(coef*x2))**4
+      log_square=exp(aa*(log(coef*x))**4+bb)
+      return
+      end function log_square
+
+      real function log_affine(x,y1,y2,x1,x2)
+      real, intent(in) :: x,y1,y2,x1,x2
+      real aa,bb
+      aa=(log(y2)-log(y1))/(log(x2)-log(x1))
+      bb=log(y2)-aa*log(x2)
+      log_affine=exp(aa*log(x)+bb)
+      return
+      end function log_affine
+
+      !! Ratio of gamma depending on T. Linear approximation
+      real function ratio_affine(temp,n)
+      implicit none
+      real, intent(in) :: temp,n
+      real aa,bb,ratio0
+      real log_affine
+      aa=(1.-0.2)/(100.-40.)
+      bb=1.-aa*100.
+      ratio0=aa*temp+bb
+      ratio_affine=log_affine(n,1.,ratio0,1e15,1e22)
+      return
+      end function ratio_affine
+
+      !! Ratio of gamma depending on T. Linear approximation
+      real function ratio_square(temp,n)
+      implicit none
+      real, intent(in) :: temp,n
+      real aa,bb,ratio0,coef
+      real log_square
+      coef=1.e-12
+      aa=(1.-0.2)/(100.-40.)
+      bb=1.-aa*100.
+      ratio0=aa*temp+bb
+      ratio_square=log_square(n,1.,ratio0,1e15,1e22,coef)
+      return
+      end function ratio_square
+
+
+
Index: trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90	(revision 3275)
@@ -257,5 +257,5 @@
 endif ! of if (nq.ge.1)
 
-!!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+!!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
 
 if (startphy_file) then
Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3275)
@@ -300,6 +300,6 @@
       real omega(ngrid,nlayer)            ! omega velocity (Pa/s, >0 when downward)
 
-      integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
-      logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
+      integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_gas, igcm_generic_ice
+      logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer
 
       real zls                       ! Solar longitude (radians).
@@ -638,5 +638,5 @@
          qsurf_hist(:,:)=qsurf(:,:)
 
-!!         call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+!!         call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
 
 !        Initialize variable for dynamical heating and zonal wind tendency diagnostic
@@ -681,5 +681,5 @@
          endif
 
-!!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+!!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
 #endif
          if (corrk) then
@@ -702,5 +702,5 @@
          endif
 
-!!         call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+!!         call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
          ! XIOS outputs
 #ifdef CPP_XIOS
@@ -711,11 +711,11 @@
 #endif
 
-!!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+!!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
 
          write(*,*) "physiq: end of firstcall"
       endif ! end of 'firstcall'
 
-!!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
-!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+!!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
+!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
       if (check_physics_inputs) then
@@ -882,6 +882,6 @@
             ! Eclipse incoming sunlight !AF24: removed
 
-!!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+!!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
+!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
 
@@ -900,12 +900,12 @@
                   do iq=1,nq
 
-                     call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-
-                     if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
+                     call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
+
+                     if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
 
                         epsi_generic=constants_epsi_generic(iq)
 
-                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
-                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
+                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_gas))
+                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_gas))
 
                      endif
@@ -931,5 +931,5 @@
                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                          &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
-                              tsurf,fract,dist_star,aerosol,muvar,                &
+                              zzlay,tsurf,fract,dist_star,aerosol,muvar,                &
                               zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
                               fluxsurfabs_sw,fluxtop_lw,                          &
@@ -1039,6 +1039,6 @@
       endif ! of if (callrad)
 
-!!      call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+!!      call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
+!!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
 
@@ -1137,5 +1137,5 @@
          ! endif
 
-!!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_gas))
 
          if (tracer) then
@@ -1144,6 +1144,6 @@
          end if ! of if (tracer)
 
-!!         call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!         call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
+!!         call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
          ! test energy conservation
@@ -1264,6 +1264,6 @@
          dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
 
-!!         call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!         call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
+!!         call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
          ! test energy conservation
@@ -1458,6 +1458,6 @@
                cloudfrac(:,:)=0.0
                do iq=1,nq
-                  call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-                  if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
+                  call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
+                  if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
                      do l = 1, nlayer
                         do ig=1,ngrid
@@ -1491,5 +1491,5 @@
             dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
 
-!!            call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!            call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
 
             ! ! Test water conservation !AF24: removed
@@ -1497,6 +1497,6 @@
          end if ! end of 'sedimentation'
 
-!!         call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!         call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
+!!         call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
   ! ---------------
@@ -1509,6 +1509,6 @@
             zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0
             !    (   zdqevap(1:ngrid,1:nlayer)                          &
-            !    !   + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
-            !    !   + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)             &
+            !    !   + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas)             &
+            !    !   + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas)             &
             !      + dqvaplscale(1:ngrid,1:nlayer) )
 
@@ -1536,8 +1536,8 @@
          endif
 
-!         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-
-!!         call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!         call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+!         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
+
+!!         call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
+!!         call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
 
@@ -1794,12 +1794,12 @@
          do iq=1,nq
 
-            call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-
-            if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
+            call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
+
+            if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
 
                do l = 1, nlayer
                   do ig=1,ngrid
                      call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq))
-                     RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_vap) / qsat_generic(ig,l,iq)
+                     RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_gas) / qsat_generic(ig,l,iq)
                   enddo
                enddo
@@ -2364,7 +2364,7 @@
 
       do iq=1,nq
-         call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-
-         if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
+         call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
+
+         if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
 
             reffrad_generic_zeros_for_wrf(:,:) = 1.
@@ -2373,5 +2373,5 @@
             comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
             comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
-            comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_vap)
+            comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas)
             comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
             ! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
@@ -2417,5 +2417,5 @@
 !        do ig=1,ngrid
 !          if ((tracer).and.(water)) then
-!           pdq(ig,:,igcm_h2o_vap) = pdq(ig,:,igcm_h2o_vap) + lsf_dq(:)
+!           pdq(ig,:,igcm_h2o_gas) = pdq(ig,:,igcm_h2o_gas) + lsf_dq(:)
 !          endif
 !          pdt(ig,:) = pdt(ig,:) + lsf_dt(:)
Index: trunk/LMDZ.PLUTO/libf/phypluto/sfluxi.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/sfluxi.F	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/sfluxi.F	(revision 3275)
@@ -1,13 +1,13 @@
       module sfluxi_mod
-      
+
       implicit none
-      
+
       contains
-      
+
       SUBROUTINE SFLUXI(PLEV,TLEV,DTAUI,TAUCUMI,UBARI,RSFI,WNOI,DWNI,
      *                  COSBI,WBARI,NFLUXTOPI,NFLUXTOPI_nu,
-     *                  FMNETI,fluxupi,fluxdni,fluxupi_nu,
+     *                  FMNETI,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,
      *                  FZEROI,TAUGSURF)
-      
+
       use radinc_h, only: NTfac, NTstart, L_LEVELS, L_NSPECTI, L_NGAUSS
       use radinc_h, only: L_NLAYRAD, L_NLEVRAD
@@ -15,12 +15,13 @@
       use comcstfi_mod, only: pi
       use gfluxi_mod, only: gfluxi
-      
+
       implicit none
-      
+
       integer NLEVRAD, L, NW, NG, NTS, NTT
-      
+
       real*8 TLEV(L_LEVELS), PLEV(L_LEVELS)
       real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
       real*8 FMNETI(L_NLAYRAD)
+      real*8 FMNETI_NU(L_NLAYRAD,L_NSPECTI)
       real*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI)
       real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
@@ -32,5 +33,5 @@
       real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)
       real*8 FTOPUP
-      
+
       real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP
       real*8 PLANCK, PLTOP
@@ -38,29 +39,30 @@
       real*8 FZEROI(L_NSPECTI)
       real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero
-      
+
       real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
       real*8 PLANCKSUM,PLANCKREF
-      
+
 ! AB : variables for interpolation
       REAL*8 C1
       REAL*8 C2
       REAL*8 P1
-      
+
 !======================================================================C
-      
+
       NLEVRAD = L_NLEVRAD
-      
+
 ! ZERO THE NET FLUXES
       NFLUXTOPI = 0.0D0
-      
+
       DO NW=1,L_NSPECTI
         NFLUXTOPI_nu(NW) = 0.0D0
         DO L=1,L_NLAYRAD
            FLUXUPI_nu(L,NW) = 0.0D0
+           FMNETI_nu(L,NW)  = 0.0
            fup_tmp(nw)=0.0D0
            fdn_tmp(nw)=0.0D0
         END DO
       END DO
-      
+
       DO L=1,L_NLAYRAD
         FMNETI(L)  = 0.0D0
@@ -68,8 +70,8 @@
         FLUXDNI(L) = 0.0D0
       END DO
-      
+
 ! WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED
 ! TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
-      
+
       TTOP  = TLEV(2)  ! JL12 why not (1) ???
       TSURF = TLEV(L_LEVELS)
@@ -86,5 +88,5 @@
       PLANCKREF = TSURF * TSURF
       PLANCKREF = sigma * PLANCKREF * PLANCKREF
-      
+
       DO NW=1,L_NSPECTI
 ! AB : PLANCKIR(NW,NTS) is replaced by P1, the linear interpolation result for a temperature TSURF
@@ -93,8 +95,8 @@
          PLANCKSUM = PLANCKSUM + P1 * DWNI(NW)
       ENDDO
-      
+
       PLANCKSUM = PLANCKREF / (PLANCKSUM * Pi)
 !JL12
-      
+
       DO 501 NW=1,L_NSPECTI
 ! SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
@@ -106,47 +108,49 @@
          BSURF = (1. - RSFI) * P1 * PLANCKSUM
          PLTOP = (1.0D0 - C2) * PLANCKIR(NW,NTT) + C2*PLANCKIR(NW,NTT+1)
-         
+
 ! If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the
 ! special Gauss point at the end.
          FZERO = FZEROI(NW)
-         
+
          IF(FZERO.ge.0.99) goto 40
-         
+
          DO NG=1,L_NGAUSS-1
-            
+
             if(TAUGSURF(NW,NG).lt. TLIMIT) then
                fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
                goto 30
             end if
-            
+
 ! SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR
 ! CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
 ! OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
-            
+
 !            TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
             TAUTOP = TAUCUMI(2,NW,NG)
             BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
-            
+
 ! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
 ! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
-! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 
-            
+! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
+
             CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
      *                TAUCUMI(1,NW,NG),
      *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
      *                BSURF,FTOPUP,FMUPI,FMDI)
-         
+
 ! NOW CALCULATE THE CUMULATIVE IR NET FLUX
             NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)
      *                * (1.0D0-FZEROI(NW))
-            
+
 ! and same thing by spectral band... (RDW)
             NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) + FTOPUP * DWNI(NW)
      *                       * GWEIGHT(NG) * (1.0D0-FZEROI(NW))
-            
+
             DO L=1,L_NLEVRAD-1
 !           CORRECT FOR THE WAVENUMBER INTERVALS
                FMNETI(L)  = FMNETI(L) + (FMUPI(L)-FMDI(L)) * DWNI(NW)
      *                    * GWEIGHT(NG)*(1.0D0-FZEROI(NW))
+               FMNETI_NU(L,NW)  = FMNETI_NU(L,NW) + (FMUPI(L)-FMDI(L))
+     *                    * DWNI(NW) * GWEIGHT(NG)*(1.0D0-FZEROI(NW))
                FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)
      *                    * (1.0D0-FZEROI(NW))
@@ -157,37 +161,39 @@
      *                          * GWEIGHT(NG) * (1.0D0 - FZEROI(NW))
             END DO
-            
+
    30       CONTINUE
-         
+
          END DO       !End NGAUSS LOOP
-         
+
    40    CONTINUE
-         
+
 ! SPECIAL 17th Gauss point
          NG     = L_NGAUSS
-         
+
 !         TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
          TAUTOP = TAUCUMI(2,NW,NG)
          BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
-         
+
 ! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
 ! CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
-! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 
-         
+! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
+
          CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG),
      *                TAUCUMI(1,NW,NG),
      *                WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP,
      *                BSURF,FTOPUP,FMUPI,FMDI)
-         
+
 ! NOW CALCULATE THE CUMULATIVE IR NET FLUX
          NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO
-         
+
 !         and same thing by spectral band... (RW)
          NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
      *      +FTOPUP*DWNI(NW)*FZERO
-         
+
          DO L=1,L_NLEVRAD-1
 ! CORRECT FOR THE WAVENUMBER INTERVALS
             FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO
+            FMNETI_nu(L,NW)  = FMNETI_nu(L,NW)+(FMUPI(L)-FMDI(L))
+     *       *DWNI(NW)*FZERO
             FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*FZERO
             FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO
@@ -196,8 +202,8 @@
      *                       + FMUPI(L) * DWNI(NW) * FZERO
          END DO
-         
+
   501 CONTINUE      !End Spectral Interval LOOP
 ! *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED****
-      
+
       END SUBROUTINE SFLUXI
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/sfluxv.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/sfluxv.F	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/sfluxv.F	(revision 3275)
@@ -2,11 +2,12 @@
 
       implicit none
-      
+
       contains
-      
+
       SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
      *                  UBAR0,STEL,NFLUXTOPV,FLUXTOPVDN,
-     *                  NFLUXOUTV_nu,NFLUXGNDV_nu,
-     *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
+     *                  NFLUXOUTV_nu,NFLUXGNDV_nu,NFLUXTOPV_nu,
+     *                  FMNETV,FMNETV_NU,FLUXUPV,FLUXDNV,
+     *                  FZEROV,taugsurf)
 
       use radinc_h, only: L_TAUMAX, L_LEVELS, L_NSPECTV, L_NGAUSS
@@ -18,4 +19,5 @@
 
       real*8 FMNETV(L_NLAYRAD)
+      real*8 FMNETV_NU(L_NLAYRAD,L_NSPECTV)
       real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
       real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
@@ -27,4 +29,5 @@
       real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
       real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN
+      real*8 NFLUXTOPV_nu(L_NSPECTV)
       real*8 NFLUXOUTV_nu(L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
@@ -50,4 +53,8 @@
          NFLUXOUTV_nu(NW)=0.0
          NFLUXGNDV_nu(NW)=0.0
+         NFLUXTOPV_nu(nw)=0.0
+         DO L=1,L_NLAYRAD
+           FMNETV_NU(L,NW) = 0.0
+        END DO
       END DO
 
@@ -64,5 +71,5 @@
 
       DO 500 NW=1,L_NSPECTV
-      
+
         F0PI = STEL(NW)
 
@@ -84,5 +91,5 @@
           BSURF = 0.
 C         LOOP OVER THE NTERMS BEGINNING HERE
- 
+
 
 !      FACTOR    = 1.0D0 - WDEL(1)*CDEL(1)**2
@@ -96,14 +103,14 @@
 C         CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
 C         WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
-C 
-C         FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 
+C
+C         FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO
 C         RETURN FLUXES FOR A GIVEN NT
 
 
           CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG),
-     *                WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW),    
+     *                WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW),
      *                BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN)
 
-C         NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 
+C         NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX
 
           NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*GWEIGHT(NG)*
@@ -114,4 +121,6 @@
             FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*
      *                           GWEIGHT(NG)*(1.0-FZEROV(NW))
+            FMNETV_NU(L,NW)=FMNETV_NU(L,NW)+( FMUPV(L)-FMDV(L) )*
+     *                           GWEIGHT(NG)*(1.0-FZEROV(NW))
             FLUXUPV(L) = FLUXUPV(L) + FMUPV(L)*GWEIGHT(NG)*
      *                   (1.0-FZEROV(NW))
@@ -120,4 +129,9 @@
           END DO
 
+c         and same thing by spectral band... (RDW)
+          NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW)
+     *      +(FLUXUP-FLUXDN)*GWEIGHT(NG)*
+     *                          (1.0-FZEROV(NW))
+
 c     band-resolved flux leaving TOA (RDW)
           NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
@@ -133,9 +147,9 @@
           DIFFVT = DIFFVT + DIFFV*GWEIGHT(NG)*(1.0-FZEROV(NW))
 
-   30     CONTINUE 
-
-        END DO   ! the Gauss loop 
-
-   40   continue 
+   30     CONTINUE
+
+        END DO   ! the Gauss loop
+
+   40   continue
 C       Special 17th Gauss point
 
@@ -143,9 +157,9 @@
 
 C       SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE VISIBLE
- 
+
         BTOP = 0.0
 
 C       LOOP OVER THE NTERMS BEGINNING HERE
- 
+
         ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX)
         BSURF = RSFV(NW)*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)
@@ -155,6 +169,6 @@
 C       CALL A SUBROUTINE THAT SOLVES  FOR THE FLUX TERMS
 C       WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER
-C 
-C       FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 
+C
+C       FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO
 C       RETURN FLUXES FOR A GIVEN NT
 
@@ -164,5 +178,5 @@
 
 
-C       NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 
+C       NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX
 
         NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*FZERO
@@ -170,8 +184,13 @@
         DO L=1,L_NLAYRAD
           FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*FZERO
+          FMNETV_NU(L,NW)=FMNETV_NU(L,NW)+( FMUPV(L)-FMDV(L) )*FZERO
           FLUXUPV(L) = FLUXUPV(L) + FMUPV(L)*FZERO
           FLUXDNV(L) = FLUXDNV(L) + FMDV(L)*FZERO
         END DO
 
+c         and same thing by spectral band... (RDW)
+          NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW)
+     *      +(FLUXUP-FLUXDN)*FZERO
+
 c     band-resolved flux leaving TOA (RDW)
           NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
@@ -196,3 +215,3 @@
 
       end module sfluxv_mod
-      
+
Index: trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90	(revision 3275)
@@ -17,16 +17,16 @@
 !$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe)
        real ttop
-       real,allocatable,dimension(:) :: kp ! TB ref pressure 
+       real,allocatable,dimension(:) :: kp ! TB ref pressure
        real p00
 !$OMP THREADPRIVATE(ttop,kp,p00)
 ! surface properties ! TB16
        real alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a
-!$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a)       
+!$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a)
        real emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq
-!$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq)       
+!$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq)
        real ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe
-!$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe)       
+!$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe)
        real alb_tho_spe,emis_tho_spe
-!$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe)       
+!$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe)
 
        real,allocatable,dimension(:) :: dryness  !"Dryness coefficient" for grnd water ice sublimation
Index: trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90	(revision 3275)
+++ trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90	(revision 3275)
@@ -0,0 +1,499 @@
+      subroutine surfprop(ngrid,nq,pfra,qsurf,tsurface,     &
+                     tidat,capcal,adjust,dist,albedo,emis,fluold, &
+                     ptimestep,zls)
+
+    !   use comgeomfi_h, only:
+      use comcstfi_mod, only: pi
+      use comsoil_h, only: nsoilmx,mlayer,volcapa,inertiedat
+      use geometry_mod, only: latitude, longitude
+    !   use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
+      use callkeys_mod, only: methane,carbox,mode_n2,mode_ch4,mode_tholins, &
+                            feedback_met,fdch4_depalb,fdch4_finalb,         &
+                            thres_ch4ice,thres_coice,thres_n2ice,           &
+                            changeti,changetid,deltab,specalb,              &
+                            fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_latn,    &
+                            fdch4_maxalb,fdch4_maxice,fdch4_ampl,           &
+                            tholateq,tholatn,tholats,tholone,tholonw,       &
+                            metlateq,metls1,metls2
+      use surfdat_h, only: alb_n2a,alb_n2b,alb_ch4,alb_ch4_eq,alb_ch4_s,    &
+                           alb_co,alb_tho,alb_tho_eq,alb_tho_spe,albedodat, &
+                           emis_n2b,emis_n2a,emis_ch4,emis_co,              &
+                           emis_tho,emis_tho_eq,emis_tho_spe,               &
+                           itch4,itch4d,ith2o,ith2od,itn2,itn2d
+      USE tracer_h, only: igcm_ch4_ice,igcm_co_ice,igcm_n2
+      use time_phylmdz_mod, only: daysec
+
+      implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     set up the properties of pluto's surface
+!
+!     Inputs
+!     ------
+!     ngrid                 Number of vertical columns
+!     nlayer                Number of layers
+!     qsurf(ngrid,iq)       Amount of ice on surface kg/m2
+!     tsurface(ngrid)       surface temperature
+!
+!     Outputs
+!     -------
+!     albedo(ngrid)
+!     emis(ngrid)
+!
+!     Both
+!     ----
+!
+!     Authors
+!     -------
+!     Tanguy Bertrand
+!
+!==================================================================
+
+!-----------------------------------------------------------------------
+!     Arguments
+
+      INTEGER ngrid, nq
+      REAL,INTENT(IN) :: qsurf(ngrid,nq)
+      REAL,INTENT(IN) :: fluold(ngrid,nq)
+      REAL,INTENT(IN) :: ptimestep
+      REAL,INTENT(IN) :: zls
+      REAL,INTENT(IN) :: tsurface(ngrid)
+      REAL,INTENT(IN) :: capcal(ngrid)
+      REAL,INTENT(IN) :: adjust
+      REAL,INTENT(IN) :: dist
+      REAL,INTENT(OUT) :: albedo(ngrid)
+      REAL,INTENT(OUT) :: emis(ngrid)
+      REAL,INTENT(OUT) :: tidat(ngrid,nsoilmx)
+      REAL,INTENT(IN) :: pfra(ngrid)
+!-----------------------------------------------------------------------
+!     Local variables
+      REAL stephan
+      DATA stephan/5.67e-08/  ! Stephan Boltzman constant
+      SAVE stephan
+      REAL tsa,tsb,coef,aa,bb
+      REAL emin,emax,emax2,tid
+      REAL n2cover,n2coversun,gamm,frac,pls,facls
+      INTEGER ig,isoil,ice_case
+      LOGICAL firstcall
+      SAVE firstcall
+      DATA firstcall/.true./
+!-----------------------------------------------------------------------
+!     Special aging for methane (not implemented)
+    !   REAL albstep            ! Time constant for albedo change (age)
+    !   DATA albstep/1.e7/
+    !   REAL albmet(ngrid)
+    !   SAVE albmet
+    !   REAL albmetmin          ! Min alb for methane ice
+    !   DATA albmetmin/0.5/
+    !   REAL albmetmax          ! Max alb for methane ice
+    !   DATA albmetmax/0.8/
+
+!-----------------------------------------------------------------------
+!  1) ALBEDOS and EMISSIVITY
+!       A. N2
+!                 CASE (0) ! fixed albedo
+!                 CASE (1) ! Albedo decreases with thickness
+!                 CASE (2) ! Special Sputnik differences of albedo
+!                 CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimationi rates
+!                 CASE (4) ! Albedo Difference in N/S (e.g. used for Triton)
+!                 CASE (5) ! Special Sputnik differences of albedo in small (1 pixel) patches (e.g. simulating dark patches / plumes)
+!            --> EMISSIVITY N2: based on the alpha/beta transition
+!       B. CO
+!       C. CH4
+!                 CASE (0) ! 2 albedos, one for the tropics, one for the poles
+!                 CASE (1) ! 3 albedos, one for the tropics, 2 for the poles (north and south)
+!                 CASE (2) ! 2 albedos + albedo feedback
+!                          SELECT CASE (feedback_met)
+!                            CASE (0) ! Default (linear from alb_ch4_eq)
+!                            CASE (1) ! Hyperbolic tangent old
+!                            CASE (2) ! hyperbolic tangent old
+!                            CASE (3) ! hyperbolic tangent equation with parameters
+!                 CASE (3) ! Eq, poles N, pole S + depending on Ls
+!       D. Tholins
+!                 CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles
+!                 CASE (1) ! Special mode one region with a different albedo
+!       E. Tholins read from file
+                 ! specalb
+
+! 2) Changes of Thermal inertia with time
+      ! if (changeti) : change of seasonal TI
+           ! if (changetid) : change of diurnal TI
+
+! 3) OTHER TESTS
+!-----------------------------------------------------------------------
+! 1) ALBEDOS and EMISSIVITY
+!-----------------------------------------------------------------------
+      ! Some input parameters:
+      pls=zls*180./pi
+      ! for equation of feedbackalbedo
+      if (methane.and.mode_ch4.eq.2.and.feedback_met.eq.3) then
+         aa=fdch4_finalb-fdch4_depalb
+         bb=fdch4_finalb
+      endif
+
+      ! Loop on all points
+      DO ig=1,ngrid
+
+        ! Looking for dominant ice:
+        ice_case=0  ! Tholins
+        if (methane) then
+          if (qsurf(ig,igcm_ch4_ice).ge.thres_ch4ice) then
+            ice_case=1
+          endif
+        endif
+        if (carbox) then
+          if (qsurf(ig,igcm_co_ice).ge.thres_coice) then
+            ice_case=2
+          endif
+        endif
+        if (qsurf(ig,igcm_n2).ge.thres_n2ice) then
+            ice_case=3
+        endif
+
+      !---------------------------------
+      ! 1.A.  N2
+      !---------------------------------
+
+        if (ice_case.eq.3) then
+
+               ! EMISSIVITY N2
+               ! change emis if phase alpha different de phase beta
+               frac=0.5*(1.+tanh(6.*(35.6-tsurface(ig))/0.5))
+               if (frac.gt.1.) frac=1.
+               emis(ig)=frac*emis_n2a+(1.-frac)*emis_n2b
+
+               ! ALBEDO N2 and emissivity changes
+               SELECT CASE (mode_n2)
+
+                 CASE (0) ! fixed albedo
+                     albedo(ig)=min(max(alb_n2b+adjust,0.),1.)
+
+                 CASE (1) ! Albedo decreases with thickness
+                     albedo(ig)=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b
+                     albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b)   ! should not be higher than alb_n2b and not lower than alb_n2b-deltab
+                 CASE (2) ! Special Sputnik differences of albedo
+                     if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)*180./pi.le.-161.or.longitude(ig)*180./pi.ge.146)) then
+                         if ( (latitude(ig)*180./pi.ge.25.).or. &
+                                ((longitude(ig)*180./pi.gt.140.).and. &
+                                 (longitude(ig)*180./pi.lt.165.)) ) then
+                                         albedo(ig)=alb_n2b-deltab
+                         else
+                                         albedo(ig)=alb_n2b+deltab
+                         endif
+                     else
+                             albedo(ig)=alb_n2b
+                     endif
+
+                 CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimation rates
+                     albedo(ig)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig)
+                     albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab)
+                     !! Triton:
+                     !albedo(ig)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig)
+                     !albedo(ig)=min(max(alb_n2b-deltab,albedo(ig)),alb_n2b+deltab)
+
+                 CASE (4) ! Albedo Difference in N/S
+                     if (latitude(ig)*180./pi.ge.0.) then
+                        albedo(ig)=min(max(alb_n2b-deltab+adjust,0.),1.)
+                     else
+                        albedo(ig)=min(max(alb_n2b+deltab+adjust,0.),1.)
+                     endif
+
+                 CASE (5) ! Special Sputnik differences of albedo in patches
+                     !! Patches Nogcm
+                     !if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)*180./pi.le.180.).and.(longitude(ig)*180./pi.ge.174.) ) then
+                     !    if (((latitude(ig)*180./pi.le.46.).and.(latitude(ig)*180./pi.ge.42.)) &
+                     ! .or.  ((latitude(ig)*180./pi.le.36.).and.(latitude(ig)*180./pi.ge.32.)) &
+                     ! .or.  ((latitude(ig)*180./pi.le.26.).and.(latitude(ig)*180./pi.ge.22.)) &
+                     ! .or.  ((latitude(ig)*180./pi.le.16.).and.(latitude(ig)*180./pi.ge.12.)) &
+                     ! .or.  ((latitude(ig)*180./pi.le.6.).and.(latitude(ig)*180./pi.ge.2.))   &
+                     !       ) then
+                     !                    albedo(ig)=alb_n2b-deltab
+
+                     !! Patches GCM
+                     if (qsurf(ig,igcm_n2).ge.1.e4) then
+                         if (((latitude(ig)*180./pi.lt.33.).and.(latitude(ig)*180./pi.gt.32.).and. &
+                              (longitude(ig)*180./pi.lt.165.5).and.(longitude(ig)*180./pi.gt.164.5)) &
+                         .or. ((latitude(ig)*180./pi.lt.30.5).and.(latitude(ig)*180./pi.gt.29.5).and. &
+                              (longitude(ig)*180./pi.lt.169.).and.(longitude(ig)*180./pi.gt.168.))) then
+                                         albedo(ig)=alb_n2b-deltab
+                         else if (((latitude(ig)*180./pi.lt.30.5).and.(latitude(ig)*180./pi.gt.29.).and. &
+                                   (longitude(ig)*180./pi.lt.165.5).and.(longitude(ig)*180./pi.gt.164.5)) &
+                             .or. ((latitude(ig)*180./pi.lt.33.).and.(latitude(ig)*180./pi.gt.32.).and.  &
+                                   (longitude(ig)*180./pi.lt.169.).and.(longitude(ig)*180./pi.gt.168.))) then
+                                         albedo(ig)=alb_n2b+deltab
+                         else
+                             albedo(ig)=alb_n2b
+                         endif
+                     else
+                             albedo(ig)=alb_n2b
+                     endif
+
+                 CASE (7) ! Albedo from albedodat and adjusted emissivity
+                   albedo(ig)=albedodat(ig)
+                   ! adjust emissivity if convergeps = true
+                   emis(ig)=min(max(emis(ig)+adjust,0.),1.)
+
+                 CASE DEFAULT
+                     write(*,*) 'STOP in surfprop.F90: mod_n2 not found'
+                     stop
+               END SELECT
+
+
+      !---------------------------------
+      ! 1.B.  CO
+      !---------------------------------
+
+        else if (ice_case.eq.2) then
+               albedo(ig)=alb_co
+               emis(ig)=emis_co
+
+      !---------------------------------
+      ! 1.C.  CH4
+      !---------------------------------
+
+        else if (ice_case.eq.1) then
+
+               SELECT CASE (mode_ch4)
+
+                 CASE (0) ! 2 albedos, one for the tropics, one for the poles
+                   emis(ig)=emis_ch4
+                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
+                      albedo(ig)=alb_ch4_eq
+                   else
+                      albedo(ig)=alb_ch4
+                   endif
+
+                 CASE (1) ! 3 albedos, one for the tropics, 2 for the poles (north and south)
+                   emis(ig)=emis_ch4
+                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
+                      albedo(ig)=alb_ch4_eq
+                   else if (latitude(ig)*180./pi.le.-metlateq) then
+                      albedo(ig)=alb_ch4_s
+                   else
+                      albedo(ig)=alb_ch4
+                   endif
+
+                 CASE (2) ! 2 albedos + albedo feedback
+
+                   emis(ig)=emis_ch4
+                   albedo(ig)=alb_ch4
+
+                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
+                      albedo(ig)=alb_ch4_eq
+                   endif
+
+                   !! Albedo feedback
+                   if (latitude(ig)*180./pi.le.fdch4_latn.and.latitude(ig)*180./pi.gt.fdch4_lats) then
+                    if (longitude(ig)*180./pi.le.fdch4_lone.and.longitude(ig)*180./pi.gt.fdch4_lonw) then
+                      if (qsurf(ig,igcm_ch4_ice).lt.fdch4_maxice) then ! fd would not apply on BTD
+                          SELECT CASE (feedback_met)
+                            CASE (0) ! Default (linear from alb_ch4_eq)
+                              albedo(ig)=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq
+                              !emis(ig)=(1.-emis_ch4)/0.002*qsurf(ig,igcm_ch4_ice)+emis_ch4
+                              if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
+                              !if (emis(ig).gt.1.) emis(ig)=1.
+                            CASE (1) ! Hyperbolic tangent old
+                              albedo(ig)=0.6*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 !
+                              !emis(ig)=0.2*0.5*(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.8 !
+                              if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
+                              !if (emis(ig).gt.1.) emis(ig)=1.
+                            CASE (2) ! hyperbolic tangent old
+                              albedo(ig)=0.5*0.6*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 !
+                              !emis(ig)=0.2*0.5*(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.8 !
+                              if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
+                              !if (emis(ig).gt.1.) emis(ig)=1.
+                            CASE (3) ! hyperbolic tangent equation with parameters
+                              albedo(ig)=aa*(-1+tanh(fdch4_ampl*qsurf(ig,igcm_ch4_ice)))+bb
+                              if (albedo(ig).gt.fdch4_maxalb) albedo(ig)=fdch4_maxalb
+                            CASE DEFAULT
+                              write(*,*) 'STOP surfprop.F90: fd wrong'
+                              stop
+                          END SELECT
+                      endif
+                    endif
+                   endif
+
+                 CASE (3) ! Eq, poles N, pole S + depending on Ls
+                   emis(ig)=emis_ch4
+                   if (latitude(ig)*180./pi.le.metlateq.and.latitude(ig)*180./pi.gt.-metlateq) then
+                      albedo(ig)=alb_ch4_eq     ! Pure methane ice
+                   else if (latitude(ig)*180./pi.le.-metlateq) then
+                      albedo(ig)=alb_ch4_s     ! Pure methane ice
+                      if (pls.le.metls2.and.pls.gt.metls1) then
+                        albedo(ig)=alb_ch4_s+deltab  ! Also used for N2, careful
+                      endif
+                   else
+                      albedo(ig)=alb_ch4     ! Pure methane ice
+                   endif
+
+                 CASE (4) ! Albedo from albedodat
+                   emis(ig)=emis_ch4
+                   albedo(ig)=albedodat(ig)
+
+                 CASE DEFAULT
+                     write(*,*) 'STOP in surfprop.F90:mod_ch4 not found'
+                     stop
+                END SELECT
+
+      !---------------------------------
+      ! 1.D.  THOLINS
+      !---------------------------------
+
+        else
+
+                SELECT CASE (mode_tholins)
+
+                 CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles
+
+                   if (latitude(ig)*180./pi.le.tholateq.and.latitude(ig)*180./pi.gt.-tholateq) then
+                      albedo(ig)=alb_tho_eq
+                      emis(ig)=emis_tho_eq
+                   else
+                      albedo(ig)=alb_tho      ! default = 0.1
+                      emis(ig)=emis_tho
+                   endif
+
+                 CASE (1) ! Special mode one region with a different albedo
+
+                   if (latitude(ig)*180./pi.le.tholateq.and.latitude(ig)*180./pi.gt.-tholateq) then
+                      albedo(ig)=alb_tho_eq
+                      emis(ig)=emis_tho_eq
+                   else
+                      albedo(ig)=alb_tho      ! default = 0.1
+                      emis(ig)=emis_tho
+                   endif
+
+                   if (latitude(ig)*180./pi.le.tholatn.and.latitude(ig)*180./pi.ge.tholats &
+                  .and.longitude(ig)*180./pi.ge.tholone.and.longitude(ig)*180./pi.le.tholonw) then
+                      albedo(ig)=alb_tho_spe
+                      emis(ig)=emis_tho_spe
+                   endif
+
+                 CASE (2) ! Depends on the albedo map
+
+                   albedo(ig)=albedodat(ig)
+                   if (albedo(ig).gt.alb_ch4) then
+                      emis(ig)=emis_ch4
+                   else
+                      emis(ig)=emis_tho
+                   endif
+
+                 CASE DEFAULT
+                     write(*,*) 'STOP in surfprop.F90:mod_ch4 not found'
+                     stop
+                END SELECT
+
+      !---------------------------------
+      ! 1.E.  Tholins read from file
+      !---------------------------------
+
+                if (specalb) then
+                  albedo(ig)=albedodat(ig)        ! Specific ground properties
+                  !if (albedodat(ig).lt.0.25) then
+                  !    albedo(ig)=alb_tho_eq
+                  !else if (albedodat(ig).lt.0.40) then
+                  !    albedo(ig)=0.35
+                  !else if (albedodat(ig).lt.0.70) then
+                  !    albedo(ig)=0.51
+                  !endif
+                endif
+
+        endif  ! ice_case
+
+      ENDDO ! ig ngrid
+
+!-----------------------------------------------------------------------
+! 2) Changes of Thermal inertia with time
+!-----------------------------------------------------------------------
+
+      IF (changeti) then ! change of seasonal TI
+        facls=8.
+        DO ig=1,ngrid
+
+           ! get depth of the different layers
+           ! limit diurnal / seasonal
+           if (changetid) then
+              if (qsurf(ig,igcm_n2)>1.e-3) then
+                 emin=facls*ITN2d/volcapa*(daysec/pi)**0.5
+                 tid=ITN2d
+              else if (qsurf(ig,igcm_ch4_ice)>1.e-3) then
+                 emin=facls*ITCH4d/volcapa*(daysec/pi)**0.5
+                 tid=ITCH4d
+              else
+                 emin=facls*ITH2Od/volcapa*(daysec/pi)**0.5
+                 tid=ITH2Od
+              endif
+           else if (ngrid.ne.1) then
+              emin=facls*inertiedat(ig,1)/volcapa*(daysec/pi)**0.5
+              tid=inertiedat(ig,1)
+           else
+              emin=facls*ITH2Od/volcapa*(daysec/pi)**0.5
+              tid=ITH2Od
+           endif
+
+           ! limit for N2
+           emax=emin+qsurf(ig,igcm_n2)/1000.
+
+           ! limit for CH4
+           if (methane) then
+               emax2=emax+qsurf(ig,igcm_ch4_ice)/1000.
+           else
+               emax2=0.
+           endif
+
+           do isoil=0,nsoilmx-1
+              if (mlayer(isoil).le.emin) then ! diurnal TI
+                   tidat(ig,isoil+1)=tid
+              else if (isoil.gt.0.and.(mlayer(isoil).gt.emin).and.(mlayer(isoil-1).lt.emin)) then ! still diurnal TI
+                   tidat(ig,isoil+1)=tid
+              else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax)) then ! TI N2
+                   tidat(ig,isoil+1)=ITN2
+              else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax2)) then
+                   tidat(ig,isoil+1)=ITCH4
+              else
+                   tidat(ig,isoil+1)=ITH2O
+              endif
+
+           enddo
+        ENDDO
+        print*, emin
+        print*, tidat(1,:)
+        print*, mlayer(:)
+
+      ELSE
+
+        DO ig=1,ngrid
+           tidat(ig,:)=inertiedat(ig,:)
+        ENDDO
+
+      ENDIF
+
+
+!-----------------------------------------------------------------------
+! 3) Tests changements emis
+!-----------------------------------------------------------------------
+           !n2cover=0.
+           !n2coversun=0.
+           !DO ig=1,ngrid
+           !   if (qsurf(ig,igcm_n2).ge.0.001) then
+           !      n2cover=n2cover+area(ig)
+           !      if (pfra(ig).gt.0.) then
+           !         n2coversun=n2coversun+area(ig)
+           !      endif
+           !    endif
+           !enddo
+           !gamm=n2cover/n2coversun
+           !gamm=1.
+           !tsb=(1/gamm*Fat1AU/dist**2*(1.-alb_n2b)/(emis_n2b*stephan))**(1/4.)
+           !tsa=(1/gamm*Fat1AU/dist**2*(1.-alb_n2b)/(emis_n2a*stephan))**(1/4.)
+           !frac=max((35.6-tsb)/abs(tsa-tsb),0.)
+           !write(*,*) 'gamm=',gamm,n2cover,n2coversun
+           !write(*,*) 'tsb=',tsb
+           !write(*,*) 'tsa=',tsa
+           !write(55,*) 'frac=',frac,tsb,tsa
+
+
+      end subroutine surfprop
+
Index: trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90	(revision 3275)
@@ -51,5 +51,5 @@
        ! Lists of constants for condensable tracers
        real, save, allocatable :: constants_mass(:)                    ! molecular mass of the specie (g/mol)
-       real, save, allocatable :: constants_delta_vapH(:)           ! Enthalpy of vaporization (J/mol)
+       real, save, allocatable :: constants_delta_gasH(:)           ! Enthalpy of vaporization (J/mol)
        real, save, allocatable :: constants_Tref(:)                 ! Ref temperature for Clausis-Clapeyron (K)
        real, save, allocatable :: constants_Pref(:)                 ! Reference pressure for Clausius Clapeyron (Pa)
@@ -58,5 +58,5 @@
        real, save, allocatable :: constants_metallicity_coeff(:)    ! Coefficient to take into account the metallicity
        real, save, allocatable :: constants_RCPV_generic(:)                   ! specific heat capacity of the tracer vapor at Tref
-!$OMP THREADPRIVATE(constants_mass,constants_delta_vapH,constants_Tref)
+!$OMP THREADPRIVATE(constants_mass,constants_delta_gasH,constants_Tref)
 !$OMP THREADPRIVATE(constants_Pref,constants_epsi_generic)
 !$OMP THREADPRIVATE(constants_RLVTT_generic,constants_metallicity_coeff,constants_RCPV_generic)
Index: trunk/LMDZ.PLUTO/libf/phypluto/vdifc_mod.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/vdifc_mod.F	(revision 3258)
+++ trunk/LMDZ.PLUTO/libf/phypluto/vdifc_mod.F	(revision 3275)
@@ -435,5 +435,5 @@
          do iq=1,nq
 
-            ! if (iq.ne.igcm_h2o_vap) then
+            ! if (iq.ne.igcm_h2o_gas) then
 
                DO ig=1,ngrid
