Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 595)
+++ trunk/LMDZ.GENERIC/README	(revision 596)
@@ -632,2 +632,9 @@
 - Added new energy conservation diagnostic variables and output
 
+== 26/03/2012 == JL
+- Added double gray case (if graybody=true in callphys.def):
+   - opacities are set to a constant value in sugas_corrk.
+   - the values are kappa_IR m^2/kg in the infrared (to be read in callphys.def)
+                    kappa_VI m^2/kg in the visible (to be read in callphys.def)
+- Cleaned continuum part in optc*
+- Added .def files for a typical 1d earth case in deftank (dry case for the moment)
Index: trunk/LMDZ.GENERIC/deftank/earth1d.gases.def
===================================================================
--- trunk/LMDZ.GENERIC/deftank/earth1d.gases.def	(revision 596)
+++ trunk/LMDZ.GENERIC/deftank/earth1d.gases.def	(revision 596)
@@ -0,0 +1,15 @@
+# gases
+4
+N2_
+CO2
+CH4
+H2O
+0.999611900E+00
+0.375999962E-03
+0.174999983E-05
+-1
+
+# First line is number of gases
+# Followed by gas names (always 3 characters)
+# and then mixing ratios.
+# mixing ratio -1 means the gas is variable.
Index: trunk/LMDZ.GENERIC/deftank/earth1d.rcm1d.def
===================================================================
--- trunk/LMDZ.GENERIC/deftank/earth1d.rcm1d.def	(revision 596)
+++ trunk/LMDZ.GENERIC/deftank/earth1d.rcm1d.def	(revision 596)
@@ -0,0 +1,102 @@
+#---------------------------------#
+# Run parameters for the 1D model #                                         
+#---------------------------------#
+
+## Planetary constants
+## ~~~~~~~~~~~~~~~~~~~~~~~
+## NB: those are mandatory
+## ~~~~~~~~~~~~~~~~~~~~~~~
+## LENGTH OF A DAY in s
+daysec         = 86400.
+## GRAVITY in m s-2
+g              = 9.81
+## LENGTH OF A YEAR in days
+year_day       = 365
+## MIN DIST STAR-PLANET in AU [periastron]
+periastr       = 1.00
+## MAX DIST STAR-PLANET in AU [apoastron]
+apoastr        = 1.00
+## DATE OF PERIASTRON in days
+peri_day       = 0.
+## OBLIQUITY in deg
+obliquit       = 0.
+## SURFACE PRESSURE in Pa
+psurf          = 101325.
+
+## Time integration parameters
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Initial date (in solar days,=0 at Ls=0)
+day0           = 0
+# Initial local time (in hours, between 0 and 24)
+time           = 12 
+# Number of time steps per sol
+day_step       = 12 
+# Number of sols to run 
+ndt            = 1300
+# Number of steps between each writting in diagfi 
+ecritphy=120
+
+
+## Vertical levels
+## ~~~~~~~~~~~~~~~
+# hybrid vertical coordinate ? (.true. for hybrid and .false. for sigma levels)
+hybrid         = .true.
+# autocompute vertical discretisation? (useful for exoplanet runs)
+autozlevs      = .false.
+# Ceiling pressure (Pa) ?
+pceil          = 0.00001
+
+## Thermal properties
+## ~~~~~~~~~~~~~~~~~~
+# Simulate global averaged conditions ? 
+global1d       = .true.
+# Latitude (deg) [only used if global1d = F]
+latitude       = 0.0
+# Solar Zenith angle (deg) [only used if global1d = T]
+szangle        = 60. 
+# Force specific heat capacity and molecular mass values
+force_cpp      = .false.
+# Specific heat capacity in J K-1 kg-1 [only used if force_cpp = T]
+cpp            = 0.
+# Molecular mass in g mol-1 [only used if force_cpp = T]
+mugaz          = 0.
+# Albedo of bare ground
+albedo         = 0.16
+# Emissivity of bare ground
+emis           = 1.0
+# Soil thermal inertia (SI)
+inertia        = 400.
+# Initial CO2 ice on the surface (kg.m-2)
+co2ice         = 0.
+
+## Wind profile
+## ~~~~~~~~~~~~
+## zonal eastward component of the geostrophic wind (m/s)
+u              = 10.
+# meridional northward component of the geostrophic wind (m/s)
+v              = 0.
+
+## Initial atmospheric temperature profile
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+# Type of initial temperature profile
+#         ichoice=1   Constant Temperature:  T=tref
+#         [Mars] ichoice=2   Savidjari profile (as Seiff but with dT/dz=cte)
+#         [Mars] ichoice=3   Lindner (polar profile)
+#         [Mars] ichoice=4   inversion
+#         [Mars] ichoice=5   Seiff  (standard profile, based on Viking entry)
+#         ichoice=6   constant T  +  gaussian perturbation (levels)
+#         ichoice=7   constant T  + gaussian perturbation (km)
+#         ichoice=8   Read in an ascii file "profile" 
+ichoice        = 1
+# Reference temperature tref (K)
+tref           = 210. 
+# Add a perturbation to profile if isin=1
+isin           = 0
+# peak of gaussian perturbation (for ichoice=6 or 7)
+pic            = 26.522
+# width of the gaussian perturbation (for ichoice=6 or 7)
+largeur        = 10
+# height of the gaussian perturbation (for ichoice=6 or 7)
+hauteur        = 30.
+
Index: trunk/LMDZ.GENERIC/deftank/earth1d.traceur.def
===================================================================
--- trunk/LMDZ.GENERIC/deftank/earth1d.traceur.def	(revision 596)
+++ trunk/LMDZ.GENERIC/deftank/earth1d.traceur.def	(revision 596)
@@ -0,0 +1,3 @@
+2
+h2o_ice
+h2o_vap
Index: trunk/LMDZ.GENERIC/deftank/earth1d_Dry.callphys.def
===================================================================
--- trunk/LMDZ.GENERIC/deftank/earth1d_Dry.callphys.def	(revision 596)
+++ trunk/LMDZ.GENERIC/deftank/earth1d_Dry.callphys.def	(revision 596)
@@ -0,0 +1,163 @@
+## Orbit / general options
+## ~~~~~~~~~~~~~~~~~~~~~~~
+# Run with or without tracer transport ?
+tracer    = .false.
+# Diurnal cycle ?  if diurnal=false, diurnally averaged solar heating
+diurnal   = .false.
+# Seasonal cycle ? if season=false, Ls stays constant, to value set in "start"
+season    = .false. 
+# Tidally resonant orbit ? must have diurnal=false, correct rotation rate in newstart
+tlocked   = .false.
+# Tidal resonance ratio ? ratio T_orbit to T_rotation
+nres      = 1
+# Write some more output on the screen ?
+lwrite    = .false.
+# Save statistics in file "stats.nc" ?
+callstats = .false.
+# Test energy conservation of model physics ?
+enertest  = .true.
+
+## Directory where external input files are
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+datadir     = /d6/jllmd/datagcm
+
+## Radiative transfer options
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call radiative transfer?
+callrad    = .true.
+# the rad. transfer is computed every "iradia" physical timestep
+iradia     = 1
+# Use blackbody for stellar spectrum ?
+stelbbody    = .false.
+# Stellar blackbody temperature ?
+stelTbb      = 2000.000
+# call multilayer correlated-k radiative transfer ?
+corrk      = .true.
+# folder in which correlated-k data is stored ?
+ corrkdir    = earth
+# corrkdir   = N2CO2poor_H2Ovar
+# corrkdir   = megaCO2
+# corrkdir   = null
+# call visible gaseous absorption in radiative transfer ?
+callgasvis = .true.
+# call continuum in radiative transfer ?
+Continuum = .false.
+# Include Rayleigh scattering in the visible ?
+rayleigh   = .true.
+# Uniform absorption coefficient in radiative transfer?
+graybody     = .false.
+# Constant absorption coefficient in visible 
+#      (in m^2/kg; only if graybody=true):
+#      tau_surf= kappa*P/g
+kappa_VI = 5.e-7
+# Constant absorption coefficient in IR
+#      (in m^2/kg; only if graybody=true):
+kappa_IR = 1.e-5 
+# Use Newtonian cooling in place of radiative transfer ?
+newtonian    = .false.
+# Radiative timescale for Newtonian cooling ? [only if newtonian = T]
+tau_relax    = 30.00000
+# Test physics timescale in 1D ?
+testradtimes = .false.
+# Characteristic planetary equilibrium (black body) temperature
+# This is used only in the aerosol radiative transfer setup. (see aerave.F)
+tplanet    = 215.
+# Output spectral OLR in 1D/3D?
+specOLR    = .true.
+# Output global radiative balance in file 'rad_bal.out' - slow for 1D!!
+meanOLR    = .false.
+# Variable gas species: Radiatively active ?
+varactive  = .false.
+# Variable gas species: Fixed vertical distribution ?
+#   (not to be used in time integration mode)
+varfixed   = .false.
+# Variable gas species: Saturation percentage value at ground ?
+satval     = .0
+# Use fixed vertical profile, 1 step, no iteration ?
+kastprof     = .false.
+# Remove lower boundary (e.g. for gas giant sims)
+noradsurf    = .false. 
+
+
+## Star type
+## ~~~~~~~~~
+startype = 1
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# The choices are:
+#
+#	startype = 1		Sol        (G2V-class main sequence)
+#	startype = 2		Ad Leo     (M-class, synthetic)
+#       startype = 3            GJ644
+#       startype = 4            HD128167
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Stellar flux at 1 AU. Examples:
+# 1366.0 W m-2		Sol today
+# 1024.5 W m-2		Sol today x 0.75 = weak early Sun
+# 18.462 W m-2		The feeble Gl581
+# 19.960 W m-2		Gl581 with e=0.38 orbital average
+Fat1AU = 1366.0
+
+## Tracer and aerosol options 
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Gravitational sedimentation of tracers (just H2O ice for now) ?
+sedimentation = .false.
+
+## Other physics options
+## ~~~~~~~~~~~~~~~~~~~~~
+# call turbulent vertical diffusion ?
+calldifv = .true.
+# use turbdiff instead of vdifc ?
+UseTurbDiff = .true.
+# call convective adjustment ?
+calladj  = .true.
+# call thermal conduction in the soil ?
+callsoil = .true.
+
+#########################################################################
+## extra non-standard definitions for Earth
+#########################################################################
+ 
+## Tracer and aerosol options 
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Fixed aerosol distributions?
+aerofixed     = .true.
+# Varying H2O cloud fraction?
+CLFvarying    = .false.
+# H2O cloud fraction?
+CLFfixval     = 0.0
+# number mixing ratio of CO2 ice particles
+Nmix_co2      = 1.e5
+# number mixing ratio of water ice particles
+Nmix_h2o      = 1.e5
+# basic dust opacity
+dusttau       = 0.0
+# water cloud pressure level (norm. by psurf)
+cloudlvl      = 0.0
+
+## Water options 
+## ~~~~~~~~~~~~~
+# Model water cycle
+water         = .false.
+# Model water cloud formation
+watercond     = .false.
+# Model water precipitation (including coagulation etc.)
+waterrain     = .false.
+# WATER: Precipitation threshold (simple scheme only) ?
+rainthreshold = 0.0011
+# Include hydrology ?
+hydrology     = .false.
+# H2O snow (and ice) albedo ?
+albedosnow    = 0.5
+# Maximum sea ice thickness ?
+maxicethick   = 0.5
+# Freezing point of seawater (degrees C) ?
+Tsaldiff      = 0.0
+# Evolve surface water sources ?
+sourceevol    = .false.
+
+## CO2 options 
+## ~~~~~~~~~~~
+# call CO2 condensation ?
+co2cond       = .false.
+# Set initial temperature profile to 1 K above CO2 condensation everywhere?
+nearco2cond   = .false.
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 595)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 596)
@@ -341,6 +341,6 @@
 
             do ig=1,ngrid
-               if(tracer)then
-               !if(tracer.and.igcm_co2_ice.gt.0)then
+               !if(tracer)then
+               if(tracer.and.igcm_co2_ice.gt.0)then
 
                   if(igcm_co2_ice.lt.1)then
Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 595)
+++ trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 596)
@@ -3,11 +3,10 @@
 ! symbols '&' in columns 73 and 6
 !
-      COMMON/callkeys/callrad,corrk,calldifv,calladj                    &
+      COMMON/callkeys/callrad,corrk,calldifv,UseTurbDiff,calladj        &
      &   , co2cond,callsoil                                             &
      &   , season,diurnal,tlocked,iradia,lwrite                         &
      &   , iaervar,iddist,topdustref,callstats,calleofdump              &
      &   , enertest                                                     &
-     &   , callgasvis                                                   &
-     &   , Continuum                                                    &
+     &   , callgasvis,Continuum, graybody                               &
      &   , Nmix_co2, Nmix_h2o                                           &
      &   , dusttau                                                      &
@@ -17,5 +16,4 @@
      &   , kastprof                                                     &
      &   , noradsurf                                                    &
-     &   , graybody                                                     &
      &   , Tstrat                                                       &
      &   , newtonian                                                    &
@@ -51,8 +49,9 @@
      &   , pceil                                                    
 
-      logical callrad,corrk,calldifv,calladj,co2cond,callsoil           &
+      logical callrad,corrk,calldifv,UseTurbDiff                        &
+     &   , calladj,co2cond,callsoil                                     &
      &   , season,diurnal,tlocked,lwrite                                &
      &   , callstats,calleofdump                                        &
-     &   , callgasvis,Continuum
+     &   , callgasvis,Continuum,graybody
 
       logical enertest
@@ -79,5 +78,4 @@
       logical CLFvarying
       logical noradsurf
-      logical graybody
 
       integer iddist
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 595)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 596)
@@ -210,4 +210,9 @@
          write(*,*) " calldifv = ",calldifv
 
+         write(*,*) "use turbdiff instead of vdifc ?"
+         UseTurbDiff=.true. ! default value
+         call getin("UseTurbDiff",UseTurbDiff)
+         write(*,*) " UseTurbDiff = ",UseTurbDiff
+
          write(*,*) "call convective adjustment ?"
          calladj=.true. ! default value
@@ -280,5 +285,5 @@
          write(*,*)" kastprof = ",kastprof
 
-         write(*,*)"Uniform absorption coefficient in IR?"
+         write(*,*)"Uniform absorption in radiative transfer?"
          graybody=.false.
          call getin("graybody",graybody)
Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 595)
+++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 596)
@@ -137,33 +137,28 @@
             DCONT = 0.0 ! continuum absorption
 
+	    if(Continuum.and.(.not.graybody))then
             ! include continua if necessary
-            wn_cont = dble(wnoi(nw))
-            T_cont  = dble(TMID(k))
-            do igas=1,ngasmx
-
-               if(gfrac(igas).eq.-1)then ! variable
-                  p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
-               else
-                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
-               endif
-
-               dtemp=0.0
-               if(gnom(igas).eq.'H2_')then
-                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
-               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
-                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
-                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-
-               endif
-
-               DCONT = DCONT + dtemp
-
-            enddo
-
-
-            DCONT = DCONT*dz(k) 
-	    
-	    if(.not.Continuum)then
-               DCONT=0.0
+               wn_cont = dble(wnoi(nw))
+               T_cont  = dble(TMID(k))
+               do igas=1,ngasmx
+
+                  if(gfrac(igas).eq.-1)then ! variable
+                     p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
+                  else
+                     p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+                  endif
+
+                  dtemp=0.0
+                  if(gnom(igas).eq.'H2_')then
+                     call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                  elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
+                     p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
+                     call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+                  endif
+
+                  DCONT = DCONT + dtemp
+               enddo
+
+               DCONT = DCONT*dz(k) 
             endif
 
@@ -213,9 +208,4 @@
                
                TAUGAS  = U(k)*ANS
-
-               if(graybody)then
-                  TAUGAS = 0.0
-                  DCONT  = U(k)*3.3e-26
-               endif
 
                TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
Index: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 595)
+++ trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 596)
@@ -132,31 +132,29 @@
             DCONT = 0.0 ! continuum absorption
 
+            if(callgasvis.and.Continuum.and.(.not.graybody))then
             ! include continua if necessary
-            wn_cont = dble(wnov(nw))
-            T_cont  = dble(TMID(k))
-            do igas=1,ngasmx
-
-               if(gfrac(igas).eq.-1)then ! variable
-                  p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
-               else
-                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
-               endif
-
-               dtemp=0.0
-               if(gnom(igas).eq.'H2_')then 
-                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
-               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
-                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
-                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-               endif
-
-               DCONT = DCONT + dtemp
-
-            enddo
-
-            DCONT = DCONT*dz(k) 
-
-            if(.not.(callgasvis.and.Continuum))then
-               DCONT=0.0
+               wn_cont = dble(wnov(nw))
+               T_cont  = dble(TMID(k))
+               do igas=1,ngasmx
+   
+                  if(gfrac(igas).eq.-1)then ! variable
+                     p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
+                  else
+                     p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+                  endif
+
+                  dtemp=0.0
+                  if(gnom(igas).eq.'H2_')then 
+                     call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                  elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
+                     p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
+                     call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+                  endif
+
+                  DCONT = DCONT + dtemp
+
+               enddo
+
+               DCONT = DCONT*dz(k) 
             endif
 
@@ -198,10 +196,4 @@
 
                TAUGAS          = U(k)*ANS
-
-               if(graybody)then
-                  TAUGAS = 0.0
-                  DCONT  = 0.0
-               endif
-
 
 
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 595)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 596)
@@ -315,5 +315,5 @@
       real dEtot, dEtots, masse, AtmToSurf_TurbFlux
       real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
-      real dEzRad(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
+      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
       real madjdE(ngridmx), lscaledE(ngridmx)
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
@@ -400,7 +400,4 @@
       logical ice_update
       save ice_update
-
-!JL12 temporary to test difference in diffusion schemes
-      logical UseTurbDiff
 
 !=======================================================================
@@ -822,4 +819,6 @@
             dEtotsSW = 0.0
             dEtotsLW = 0.0
+	    dEzRadsw(:,:)=0.
+	    dEzRadlw(:,:)=0.
             do ig = 1, ngrid
                do l = 1, nlayer
@@ -827,4 +826,6 @@
                   dEtotSW = dEtotSW + cpp*masse*zdtsw(ig,l)*area(ig)
                   dEtotLW = dEtotLW + cpp*masse*zdtlw(ig,l)*area(ig)
+		  dEzRadsw(ig,l)=cpp*masse*zdtsw(ig,l)
+		  dEzRadlw(ig,l)=cpp*masse*zdtlw(ig,l)
                enddo
                dEtotsSW = dEtotsSW + fluxsurf_sw(ig)*(1.-albedo(ig))*area(ig)
@@ -862,5 +863,4 @@
 
 !JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff	 
-	 UseTurbDiff=.true.
 	 if (UseTurbDiff) then
 	 
@@ -2008,5 +2008,9 @@
 	
 	if(enertest) then
-	  call writediagfi(ngridmx,"dEzdif","atm vdifc energy change","w.m^-2",3,dEzdif)
+	  if (calldifv) call writediagfi(ngridmx,"dEzdif","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdif)
+	  if (corrk) then
+	     call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
+	     call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
+	  endif 
           if(watercond) then
 	     call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
Index: trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 595)
+++ trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 596)
@@ -12,4 +12,5 @@
 !     -------
 !     Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009)
+!     Added double gray case by Jeremy Leconte (2012)
 !
 !     Summary
@@ -25,8 +26,9 @@
       use datafile_mod, only: datadir
       use gases_h
+      use ioipsl_getincom 
       implicit none
 
 #include "callkeys.h"
-
+#include "comcstfi.h"
 !==================================================================
 
@@ -44,4 +46,5 @@
 
       real*8 x, xi(4), yi(4), ans, p
+      real kappa_IR, kappa_VI
 
       integer ngas, igas
@@ -53,5 +56,4 @@
       file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
       file_path=TRIM(datadir)//TRIM(file_id)
-
 
       ! check that the file exists
@@ -266,11 +268,9 @@
 !-----------------------------------------------------------------------
 
-
-
 !=======================================================================
 !     Get gaseous k-coefficients and interpolate onto finer pressure grid
 
       ! VISIBLE
-      if (callgasvis.and..not.TRIM(corrkdir).eq.'null') then
+      if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 
          file_path=TRIM(datadir)//TRIM(file_id)
@@ -288,5 +288,13 @@
          read(111,*) gasv8
          close(111)
-         
+	 
+      else if (callgasvis.and.graybody) then
+!    constant absorption coefficient in visible
+         write(*,*)"constant absorption coefficient in visible:"
+         kappa_VI=-100000.
+         call getin("kappa_VI",kappa_VI)
+         write(*,*)" kappa_VI = ",kappa_VI
+	 kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule         
+         gasv8(:,:,:,:,:)=kappa_VI
       else
          print*,'Visible corrk gaseous absorption is set to zero.'
@@ -295,5 +303,5 @@
 
       ! INFRA-RED
-      if (.not.TRIM(corrkdir).eq.'null') then
+      if ((.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 
          file_path=TRIM(datadir)//TRIM(file_id)
@@ -311,46 +319,26 @@
          read(111,*) gasi8
          close(111)
-     
-         do nw=1,L_NSPECTI
-            fzeroi(nw) = 0.
-!            do nt=1,L_NTREF
-!               do np=1,L_NPREF
-!                  do nh=1,L_REFVAR
-!                     do ng = 1,L_NGAUSS
-!                        if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzeroi(nw)=fzeroi(nw)+1.
-!                        endif
-!                     end do
-!                  end do
-!               end do
-!            end do
-!            fzeroi(nw)=fzeroi(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
-         end do
-
-         do nw=1,L_NSPECTV
-            fzerov(nw) = 0.
-!            do nt=1,L_NTREF
-!               do np=1,L_NPREF
-!                  do nh=1,L_REFVAR
-!                     do ng = 1,L_NGAUSS
-!                        if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzerov(nw)=fzerov(nw)+1.
-!                        endif
-!                     end do
-!                  end do
-!               end do
-!            end do
-!            fzerov(nw)=fzerov(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
-         end do
-
-
-         do nw=1,L_NSPECTV
-            fzerov(nw) = 0.
-         end do
-
+
+      else if (graybody) then
+!     constant absorption coefficient in IR
+         write(*,*)"constant absorption coefficient in InfraRed:"
+         kappa_IR=-100000.
+         call getin("kappa_IR",kappa_IR)
+         write(*,*)" kappa_IR = ",kappa_IR	 
+         kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule       
+         gasi8(:,:,:,:,:)=kappa_IR
       else
          print*,'Infrared corrk gaseous absorption is set to zero.'
          gasi8(:,:,:,:,:)=0.0
       endif
+
+      do nw=1,L_NSPECTI
+         fzeroi(nw) = 0.
+      end do
+
+      do nw=1,L_NSPECTV
+         fzerov(nw) = 0.
+      end do
+
 
 !     Take log10 of the values - this is what we will interpolate.
