Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 365)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 366)
@@ -35,4 +35,5 @@
 #include "callkeys.h"
 #include "tracer.h"
+#include "gases.h"
 
 !-----------------------------------------------------------------------
@@ -104,6 +105,6 @@
       REAL*8 tauaero(L_LEVELS+1,naerkind)
       REAL*8 nfluxtopv,nfluxtopi,nfluxtop
-      real*8 NFLUXTOPV_nu(L_NSPECTV)
-      real*8 NFLUXTOPI_nu(L_NSPECTI)
+      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)
@@ -163,5 +164,6 @@
       logical OLRz
       real OLR_nu(ngrid,L_NSPECTI)
-      real GSR_nu(ngrid,L_NSPECTV)
+      !real GSR_nu(ngrid,L_NSPECTV)
+      real OSR_nu(ngrid,L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
 
@@ -184,5 +186,5 @@
       external CBRT
 
-!     included by RW for runway greenhouse 1D study
+!     included by RW for runaway greenhouse 1D study
       real muvar(ngridmx,nlayermx+1)
       real vtmp(nlayermx)
@@ -589,5 +591,23 @@
       end do
 
+
+
+!-----------------------------------------------------------------------
+!     kcm mode only
       if(kastprof)then
+
+         DO l=1,nlayer
+            muvarrad(2*l)   = mugaz
+            muvarrad(2*l+1) = mugaz
+         END DO
+
+         do k=1,L_LEVELS
+            qvar(k) = 0.0
+         end do
+         print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
+      endif
+
+
+      if(kastprof.and.(ngasmx.gt.1))then
 
          DO l=1,nlayer
@@ -613,10 +633,7 @@
          qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
 
-         !do k=1,L_LEVELS
-         !   qvar(k) = 1.0
-         !end do
-         !print*,'ASSUMING qH2O=1 EVERYWHERE IN CALLCORRK!!'
-
-      endif
+
+      endif
+
 
       ! Keep values inside limits for which we have radiative transfer coefficients
@@ -720,9 +737,13 @@
 
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
-                 acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu,  &
+                 acosz,stel_fract,gweight,                         &
+                 nfluxtopv,nfluxoutv_nu,nfluxgndv_nu,              &
+                 !acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu,  &
                  fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
 
          else                          ! during the night, fluxes = 0
-            nfluxtopv=0.0
+            nfluxtopv       = 0.0
+            nfluxoutv_nu(:) = 0.0
+            nfluxgndv_nu(:) = 0.0
             do l=1,L_NLAYRAD
                fmnetv(l)=0.0
@@ -775,5 +796,6 @@
             end do
             do nw=1,L_NSPECTV 
-               GSR_nu(ig,nw)=nfluxgndv_nu(nw)
+               !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
+               OSR_nu(ig,nw)=nfluxoutv_nu(nw)
             end do
          endif
@@ -824,6 +846,6 @@
       if(specOLR)then
         if(ngrid.ne.1)then
-          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
-          call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W m^-2",3,GSR_nu)
+          !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
+          !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
         endif
       endif
@@ -845,7 +867,7 @@
                close(117)
 
-               open(127,file='GSRnu.out')
+               open(127,file='OSRnu.out')
                do nw=1,L_NSPECTV
-                  write(127,*) GSR_nu(1,nw)
+                  write(127,*) OSR_nu(1,nw)
                enddo
                close(127)
Index: /trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90	(revision 365)
+++ /trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90	(revision 366)
@@ -42,5 +42,10 @@
      Pn = rho_n*T*rmn
 
-     if(gnom(2).eq.'H2O')then
+     if(ngasmx.eq.1)then
+        print*,'Cannot have moist adiabat with one gas...'
+        stop
+     endif
+
+     if(gnom(ngasmx).eq.'H2O')then
 
         call psat_H2O(T-2d-1,psat_minus)
@@ -68,5 +73,5 @@
         endif
 
-     elseif(gnom(2).eq.'NH3')then
+     elseif(gnom(ngasmx).eq.'NH3')then
 
         call psat_NH3(T-2d-1,psat_minus)
@@ -104,9 +109,10 @@
   case(0) ! dry
 
-     if(gnom(2).eq.'H2O')then
+     cp_v=0.0
+     if(gnom(ngasmx).eq.'H2O')then
         cp_v = (32.24+1.923d-3*T+1.055d-5*T**2-3.511d-9*T**3)/m_v
-     elseif(gnom(2).eq.'NH3')then
+     elseif(gnom(ngasmx).eq.'NH3')then
         cp_v = 2.058d3
-     elseif(gnom(2).eq.'CH4')then
+     elseif(gnom(ngasmx).eq.'CH4')then
         cp_v = 2.226d3
      endif
Index: /trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 365)
+++ /trunk/LMDZ.GENERIC/libf/phystd/kcm1d.F90	(revision 366)
@@ -162,4 +162,5 @@
   call su_gases
   call calc_cpp_mugaz
+
 
   call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp)
@@ -192,12 +193,15 @@
            read(90,*,iostat=ierr) tnom(iq)
            if (ierr.ne.0) then
-              write(*,*) 'rcm1d: error reading tracer names...'
+              write(*,*) 'kcm1d: error reading tracer names...'
               stop
            endif
         enddo !of do iq=1,nq
      endif
-  endif
+  !endif
   
   call initracer()
+
+  endif
+
 
   do iq=1,nqmx
@@ -215,5 +219,5 @@
 
   iter    = 1
-  Tstrat  = 200.0
+  Tstrat  = 60.0
   dTstrat = 1000.0
 
@@ -227,7 +231,4 @@
      call kcmprof_fn(psurf,qsurf(1),tsurf,    &
           tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
-
-
-     !if(1.eq.2)then
 
 !    Run radiative transfer
@@ -242,10 +243,9 @@
      print*,'Tstrat = ',Tstrat
      dTstrat = Tstrat
-     Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
+     !Tstrat  = Tsurf*(0.2786*(psurf/100000.)**(-1.123))**0.25 
+     ! skin temperature (gray approx.) using analytic pure H2 expression
+     !Tstrat  = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.)
+     Tstrat  = (fluxtop_lw/(2*sigma))**0.25 ! skin temperature (gray approx.)
      dTstrat = dTstrat-Tstrat
-
-     !endif
-
-     !dTstrat = 0
 
      if(abs(dTstrat).lt.1.0)then
@@ -264,6 +264,7 @@
 
 ! Calculate total atmospheric energy
-  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
-     q(:,1),muvar,Eatmtot)
+  Eatmtot=0.0
+!  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
+!     q(:,1),muvar,Eatmtot)
 
 ! ------------------------
Index: /trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 365)
+++ /trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90	(revision 366)
@@ -41,5 +41,6 @@
   double precision Ptop, dlogp, Psat_max
   parameter (Ptop=1.0)           ! Pressure at TOA [Pa]
-  parameter (Psat_max=5000.0) ! Maximum vapour pressure [Pa]
+  !parameter (Psat_max=100000.0) ! Maximum vapour pressure [Pa]
+  parameter (Psat_max=0.0) ! set to zero for dry calculations
 
   double precision T(1:nlay)                       ! temperature [K]
@@ -83,11 +84,22 @@
   ! modify/generalise later??
 
-  if(gnom(2).eq.'H2O')then
+  if(ngasmx.gt.2)then
+     print*,'Only two species possible at present'
+     stop
+  elseif(ngasmx.eq.1)then
+     if(psat_max.gt.0.0)then
+        print*,'Must have Psat_max=0 if no variable species'
+        stop
+     endif
+     print*, 'Assuming pure atmosphere'
+     m_v   = 1.0
+     tcrit = 1000.0
+  elseif(gnom(ngasmx).eq.'H2O')then
      m_v   = dble(mH2O/1000.)
      tcrit = 6.47d2
-  elseif(gnom(2).eq.'NH3')then
+  elseif(gnom(ngasmx).eq.'NH3')then
      m_v   = 17.031/1000.
      tcrit = 4.06d2
-  elseif(gnom(2).eq.'CH4')then
+  elseif(gnom(ngasmx).eq.'CH4')then
      m_v   = 16.04/1000.
      tcrit = 1.91d2
@@ -113,11 +125,13 @@
 !i!  endif
 
-  if(gnom(2).eq.'H2O')then
+
+  psat_v=psat_max
+  if(gnom(ngasmx).eq.'H2O')then
      call Psat_H2O(tsurf,psat_v)
-  elseif(gnom(2).eq.'NH3')then
+  elseif(gnom(ngasmx).eq.'NH3')then
      call Psat_NH3(tsurf,psat_v)
   endif
 
-  ! Moist adiabat unless greater than psat_max
+  ! Moist adiabat unless greater than or equal to psat_max
   if(psat_v*1d6.lt.psat_max)then
     Psurf_v = Psat_v*1d6
@@ -141,4 +155,6 @@
   endif
 
+
+
   ! define fine pressure grid
   psurf_rcm = real(Psurf_n+Psurf_v)
@@ -211,7 +227,8 @@
 
      ! test for moist adiabat at next level
-     if(gnom(2).eq.'H2O')then
+     psat_v=psat_max
+     if(gnom(ngasmx).eq.'H2O')then
         call Psat_H2O(T(ilay+1),psat_v)
-     elseif(gnom(2).eq.'NH3')then
+     elseif(gnom(ngasmx).eq.'NH3')then
         call Psat_NH3(T(ilay+1),psat_v)
      endif
@@ -234,8 +251,9 @@
 
      if(profil_flag(ilay) .eq. 1)then
-        
-        if(gnom(2).eq.'H2O')then
+       
+        psat_v=psat_max 
+        if(gnom(ngasmx).eq.'H2O')then
            call Psat_H2O(T(ilay+1),psat_v)
-        elseif(gnom(2).eq.'NH3')then
+        elseif(gnom(ngasmx).eq.'NH3')then
            call Psat_NH3(T(ilay+1),psat_v)
         endif
Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 365)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 366)
@@ -157,7 +157,7 @@
             DCONT = DCONT*dz(k) 
 
-               if((.not.callgasvis))then
-                  DCONT=0.0
-               endif
+            if((.not.callgasvis))then
+               DCONT=0.0
+            endif
 
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F	(revision 365)
+++ /trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F	(revision 366)
@@ -1,4 +1,5 @@
       SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
-     *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUXGNDV_nu,
+     *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUXOUTV_nu,
+     *                  NFLUXGNDV_nu,
      *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
 
@@ -18,5 +19,5 @@
       real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
       real*8 NFLUXTOPV, FLUXUP, FLUXDN
-      real*8 NFLUXTOPV_nu(L_NSPECTV)
+      real*8 NFLUXOUTV_nu(L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
       real*8 GWEIGHT(L_NGAUSS)
@@ -38,5 +39,5 @@
 
       DO NW=1,L_NSPECTV
-         NFLUXTOPV_nu(NW)=0.0
+         NFLUXOUTV_nu(NW)=0.0
          NFLUXGNDV_nu(NW)=0.0
       END DO
@@ -108,11 +109,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)
+     *      +FLUXUP*GWEIGHT(NG)*(1.0-FZEROV(NW))
 
-
-c     flux at gnd (RDW)
+c     band-resolved flux at ground (RDW)
           NFLUXGNDV_nu(NW) = NFLUXGNDV_nu(NW)
      *      +FMDV(L_NLAYRAD)*GWEIGHT(NG)*(1.0-FZEROV(NW))
@@ -163,9 +162,9 @@
         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)
+     *      +FLUXUP*FZERO
 
-c     flux at gnd (RDW)
+c     band-resolved flux at ground (RDW)
           NFLUXGNDV_nu(NW) = NFLUXGNDV_nu(NW)+FMDV(L_NLAYRAD)*FZERO
 
