Index: /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 526)
@@ -235,7 +235,7 @@
                           (  0.75 * QREFvis3d(ig,l,iaer) /        &
                           ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
-                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
-                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   &
-                            CLF1
+                            pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
+                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the 
+                            CLF1				    !   clear=false/pq=0 case
 
                      aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0)) 
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 526)
@@ -1,8 +1,9 @@
-      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,    &
+      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
           albedo,emis,mu0,pplev,pplay,pt,                      & 
           tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
-          reffrad,tau_col,cloudfrac,totcloudfrac,   &
+	  OLR_nu,OSR_nu,                                       &
+          reffrad,tau_col,cloudfrac,totcloudfrac,              &
           clearsky,firstcall,lastcall)
 
@@ -83,5 +84,6 @@
       REAL fluxabs_sw(ngridmx)    ! SW flux absorbed by planet (W/m2)
       REAL fluxtop_dn(ngridmx)    ! incident top of atmosphere SW flux (W/m2)
-
+      REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
+      REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
 !-----------------------------------------------------------------------
 !     Declaration of the variables required by correlated-k subroutines
@@ -154,16 +156,6 @@
 !      real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
 
-!     for OLR spec
-      integer OLRcount
-      save OLRcount
-      integer OLRcount2
-      save OLRcount2
-      character(len=2)  :: tempOLR
-      character(len=30) :: filenomOLR
       !real ptime, pday
       logical OLRz
-      real OLR_nu(ngrid,L_NSPECTI)
-      !real GSR_nu(ngrid,L_NSPECTV)
-      real OSR_nu(ngrid,L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
 
@@ -272,4 +264,7 @@
             stop
          endif
+	 
+	 OLR_nu=0.
+	 OSR_nu=0.
 
          firstcall=.false.   
@@ -331,5 +326,5 @@
 
             do ig=1,ngrid
-               if(tracer)then
+               if(tracer.and.igcm_co2_ice.gt.0)then
                   reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
                        (4*Nmix_co2*pi*rho_co2) )
@@ -561,5 +556,5 @@
          end do
          qvar(1)=qvar(2)
-         qvar(2*nlayermx+1)=qsurf(ig,i_var)
+!         qvar(2*nlayermx+1)=qsurf(ig,i_var) !JL12 not very good to compare kg/kg and kg/m2???
 
       elseif(varfixed)then
@@ -590,7 +585,12 @@
          Ttemp = tsurf(ig)
          call watersat(Ttemp,ptemp,qsat)
+	 
 
          !qvar(2*nlayermx+1)=qsat      ! fully saturated everywhere
          qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!  JL: for completely constant water profile uncoment the following line
+!	 qvar=0.005
 
       else
@@ -766,5 +766,5 @@
               qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
               taugsurfi,qvar,muvarrad)
-
+	      
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
               wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, & 
@@ -797,9 +797,9 @@
          if(specOLR)then
             do nw=1,L_NSPECTI 
-               OLR_nu(ig,nw)=nfluxtopi_nu(nw)
+               OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
             end do
             do nw=1,L_NSPECTV 
                !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
-               OSR_nu(ig,nw)=nfluxoutv_nu(nw)
+               OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
             end do
          endif
@@ -848,12 +848,7 @@
 
 !     IR spectral output, for exoplanet observational comparison
-      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,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
-        endif
-      endif
-
-      if(lastcall.and.(ngrid.eq.1))then
+
+
+      if(lastcall.and.(ngrid.eq.1))then  ! could disable the 1D output, they are in the diagfi and diagspec... JL12
 
            print*,'Saving scalar quantities in surf_vals.out...'
@@ -864,17 +859,17 @@
            close(116)
 
-           if(specOLR)then
-               open(117,file='OLRnu.out')
-               do nw=1,L_NSPECTI
-                  write(117,*) OLR_nu(1,nw)
-               enddo
-               close(117)
-
-               open(127,file='OSRnu.out')
-               do nw=1,L_NSPECTV
-                  write(127,*) OSR_nu(1,nw)
-               enddo
-               close(127)
-           endif
+!           if(specOLR)then
+!               open(117,file='OLRnu.out')
+!               do nw=1,L_NSPECTI
+!                  write(117,*) OLR_nu(1,nw)
+!               enddo
+!               close(117)
+!
+!               open(127,file='OSRnu.out')
+!               do nw=1,L_NSPECTV
+!                  write(127,*) OSR_nu(1,nw)
+!               enddo
+!               close(127)
+!           endif
 
 !     OLR vs altitude: do it as a .txt file
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 526)
@@ -9,4 +9,5 @@
      &   , enertest                                                     &
      &   , callgasvis                                                   &
+     &   , Continuum                                                    &
      &   , Nmix_co2, Nmix_h2o                                           &
      &   , dusttau                                                      &
@@ -52,5 +53,5 @@
      &   , season,diurnal,tlocked,lwrite                                &
      &   , callstats,calleofdump                                        &
-     &   , callgasvis
+     &   , callgasvis,Continuum
 
       logical enertest
Index: /trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 526)
@@ -46,4 +46,5 @@
 !     Francois Forget (1996)
 !     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
+!     Includes simplifed nucleation by J. Leconte (2011)
 !     
 !==================================================================
@@ -96,4 +97,5 @@
       REAL zcpi
       REAL ztcond (ngridmx,nlayermx)
+      REAL ztnuc (ngridmx,nlayermx)
       REAL ztcondsol(ngridmx) 
       REAL zdiceco2(ngridmx)
@@ -271,9 +273,10 @@
 !      end if
 
-!     Forecast the atmospheric frost temperature 'ztcond'
+!     Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'
       DO l=1,nlayer
          DO ig=1,ngrid
             ppco2=gfrac(igasco2)*pplay(ig,l)
             call get_tcond_co2(ppco2,ztcond(ig,l))
+            call get_tnuc_co2(ppco2,ztnuc(ig,l))
          ENDDO
       ENDDO
@@ -379,5 +382,6 @@
 
 
-               IF ((zt(ig,l).LT.ztcond(ig,l)).or.(zq(ig,l,i_co2ice).gt.0)) THEN
+   !     ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
+               IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN 
                   condsub(ig)=.true.
                   pdtc(ig,l)   = (ztcond(ig,l) - zt(ig,l))/subptimestep
@@ -524,5 +528,5 @@
       real, parameter :: ptriple=518000.0
 
-      peff=p!/co2supsat
+      peff=p
 
       if(peff.lt.ptriple)then
@@ -535,2 +539,29 @@
 
     end subroutine get_tcond_co2
+    
+    
+    
+ 
+!-------------------------------------------------------------------------
+    subroutine get_tnuc_co2(p,tnuc)
+!   Calculates the nucleation temperature for CO2, based on a simple super saturation criterion 
+!     (JL 2011)
+
+      implicit none
+
+#include "callkeys.h"
+
+      real p, peff, tnuc
+      real, parameter :: ptriple=518000.0
+
+      peff=p/co2supsat
+
+      if(peff.lt.ptriple)then
+         tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
+      else
+         tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2 
+         ! liquid-vapour transition (based on CRC handbook 2003 data)
+      endif
+      return
+
+    end subroutine get_tnuc_co2
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 526)
@@ -192,4 +192,10 @@
          call getin("callgasvis",callgasvis)
          write(*,*) " callgasvis = ",callgasvis
+
+         write(*,*) "call continuum opacities in radiative transfer ?",
+     &              "(matters only if callrad=T)"
+         Continuum=.false. ! default value
+         call getin("Continuum",Continuum)
+         write(*,*) " Continuum = ",Continuum
          
          write(*,*) "call turbulent vertical diffusion ?"
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 526)
@@ -2,5 +2,5 @@
 
       use radinc_h, only: L_NSPECTI
-      use radcommon_h, only: WAVEI
+      use radcommon_h, only: WNOI
 
       implicit none
@@ -54,5 +54,5 @@
       integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
  !     integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
-      integer :: idim_bandsIR ! "bandsIR" dimension ID #
+      integer :: idim_bandsIR ! "IR Wavenumber" dimension ID #
       integer, dimension(2) :: id  
 
@@ -114,5 +114,5 @@
       ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
 !      ierr = NF_DEF_DIM (nid, "subsurface_layers",nsoilmx,idim_nsoilmx)
-      ierr = NF_DEF_DIM (nid, "bandsIR",L_NSPECTI,idim_bandsIR)
+      ierr = NF_DEF_DIM (nid, "IR Wavenumber",L_NSPECTI,idim_bandsIR)
 
       ierr = NF_ENDDEF(nid)
@@ -290,17 +290,19 @@
       ! define variable
 #ifdef NC_DOUBLE
-      ierr=NF_DEF_VAR(nid,"bandsIR",NF_DOUBLE,1,idim_bandsIR,nvarid)
-#else
-      ierr=NF_DEF_VAR(nid,"bandsIR",NF_FLOAT,1,idim_bandsIR,nvarid)
-#endif
-      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 27,
-     .                        "Band limits in the infrared")
-      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"cm^-1")
+      ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_DOUBLE,1,
+     .				idim_bandsIR,nvarid)
+#else
+      ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_FLOAT,1,
+     .				idim_bandsIR,nvarid)
+#endif
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 34,
+     .                        "Band mid frequency in the infrared")
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",5,"cm^-1")
       ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
       ! write variable
 #ifdef NC_DOUBLE
-      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WAVEI))
-#else
-      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WAVEI))
+      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WNOI))
+#else
+      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WNOI))
 #endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 526)
@@ -2,5 +2,5 @@
 
       use radinc_h, only: L_NSPECTV
-      use radcommon_h, only: WAVEV
+      use radcommon_h, only: WNOV
 
       implicit none
@@ -54,5 +54,5 @@
       integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
  !     integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
-      integer :: idim_bandsVI ! "bandsIR" dimension ID #
+      integer :: idim_bandsVI ! "VI Wavenumber" dimension ID #
       integer, dimension(2) :: id  
 
@@ -114,5 +114,5 @@
       ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
 !      ierr = NF_DEF_DIM (nid, "subsurface_layers",nsoilmx,idim_nsoilmx)
-      ierr = NF_DEF_DIM (nid, "bandsVI",L_NSPECTV,idim_bandsVI)
+      ierr = NF_DEF_DIM (nid, "VI Wavenumber",L_NSPECTV,idim_bandsVI)
 
       ierr = NF_ENDDEF(nid)
@@ -290,17 +290,19 @@
       ! define variable
 #ifdef NC_DOUBLE
-      ierr=NF_DEF_VAR(nid,"bandsVI",NF_DOUBLE,1,idim_bandsVI,nvarid)
-#else
-      ierr=NF_DEF_VAR(nid,"bandsVI",NF_FLOAT,1,idim_bandsVI,nvarid)
-#endif
-      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 27,
-     .                        "Band limits in the visible")
-      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"cm^-1")
+      ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_DOUBLE,1,
+     .				idim_bandsVI,nvarid)
+#else
+      ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_FLOAT,1,
+     .				idim_bandsVI,nvarid)
+#endif
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 33,
+     .                        "Band mid frequency in the visible")
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",5,"cm^-1")
       ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
       ! write variable
 #ifdef NC_DOUBLE
-      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WAVEV))
-#else
-      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WAVEV))
+      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,dble(WNOV))
+#else
+      ierr=NF_PUT_VAR_REAL (nid,nvarid,real(WNOV))
 #endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90	(revision 526)
@@ -150,5 +150,5 @@
 
       else
-  
+
          call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
          call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
@@ -162,4 +162,5 @@
          abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
          abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account
+
 
          ! unlike for Rayleigh scattering, we do not currently weight by the BB function
Index: /trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 526)
@@ -163,4 +163,8 @@
 
             DCONT = DCONT*dz(k) 
+	    
+	    if(.not.Continuum)then
+               DCONT=0.0
+            endif
 
             !--- Kasting's CIA ----------------------------------------
@@ -217,5 +221,4 @@
                TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
                !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
-
                DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 526)
@@ -157,8 +157,7 @@
             DCONT = DCONT*dz(k) 
 
-            if((.not.callgasvis))then
+            if(.not.(callgasvis.and.Continuum))then
                DCONT=0.0
             endif
-
 
             do NG=1,L_NGAUSS-1
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 526)
@@ -7,5 +7,5 @@
                   pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
  
-      use radinc_h, only : naerkind
+      use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV
       use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
       use gases_h
@@ -116,5 +116,4 @@
 #include "comsaison.h"
 #include "control.h"
-#include "comg1d.h"
 #include "tracer.h"
 #include "watercap.h"
@@ -190,5 +189,5 @@
 
       character*80 fichier 
-      integer l,ig,ierr,iq,i, tapphys
+      integer l,ig,ierr,iq,i, tapphys,nw
 
       real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
@@ -199,5 +198,7 @@
       real fluxtop_dn(ngridmx)
       real fluxdyn(ngridmx)          ! horizontal heat transport by dynamics      
-      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn
+      real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
+      real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
+      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
 
 
@@ -326,7 +327,4 @@
       real zpopskNI(ngridmx,nlayermx)
 
-!     included by RW to make 1D saves not every timestep
-      integer countG1D
-      save countG1D
 
 !     included by BC for evaporation     
@@ -357,12 +355,13 @@
 !     included by BC for double radiative transfer call
       logical clearsky
-
-      real zdtsw1(ngridmx,nlayermx),zdtsw2(ngridmx,nlayermx)
-      real zdtlw1(ngridmx,nlayermx),zdtlw2(ngridmx,nlayermx)
-      real fluxsurf_lw1(ngridmx), fluxsurf_lw2(ngridmx)     
-      real fluxsurf_sw1(ngridmx), fluxsurf_sw2(ngridmx)    
-      real fluxtop_lw1(ngridmx), fluxtop_lw2(ngridmx)     
-      real fluxabs_sw1(ngridmx), fluxabs_sw2(ngridmx)   
-      real tau_col1(ngrid), tau_col2(ngrid)
+      real zdtsw1(ngridmx,nlayermx)
+      real zdtlw1(ngridmx,nlayermx)
+      real fluxsurf_lw1(ngridmx)     
+      real fluxsurf_sw1(ngridmx)  
+      real fluxtop_lw1(ngridmx)    
+      real fluxabs_sw1(ngridmx)  
+      real tau_col1(ngrid)
+      real OLR_nu1(ngrid,L_NSPECTI)
+      real OSR_nu1(ngrid,L_NSPECTV)
       real tf, ntf
 
@@ -423,10 +422,4 @@
       if (firstcall) then
 
-         if(ngrid.eq.1)then
-            saveG1D=day_step*10
-            !saveG1D=1
-            !saveG1D=-1 ! to never save 
-            countG1D=1
-         endif
 
 !        variables set to 0
@@ -447,5 +440,5 @@
 !        read startfi (initial parameters)
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+	  print*,'in physiqu i am in phyetat0'
          call phyetat0("startfi.nc",0,0,nsoilmx,nq,           &
                day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
@@ -681,5 +674,5 @@
 
       if (callrad) then
-         if( mod(icount-1,iradia).eq.0) then
+         if( mod(icount-1,iradia).eq.0.or.lastcall) then
 
 !          Compute local stellar zenith angles
@@ -720,12 +713,22 @@
               endif
               muvar(:,:)=0.0 ! only used for climate evolution for now
-
-!             Option to call scheme once for clear regions, once for cloudy regions (BC)
+	      
+	      
+!             standard callcorrk
+              clearsky=.false.
+              call callcorrk(ngrid,nlayer,pq,nq,qsurf,         &
+                      albedo,emis,mu0,pplev,pplay,pt,                    &
+                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
+                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
+		      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
+                      reffrad,tau_col,cloudfrac,totcloudfrac,            &
+                      clearsky,firstcall,lastcall)      
+	      
+
+!             Option to call scheme once more for clear regions
               if(CLFvarying)then
 
-               !!! THIS IS NOT A VERY GOOD STRATEGY HERE (everything should be in callcorrk)
-               !!! ---> PROBLEMS WITH ALLOCATED ARRAYS 
-               !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
-
+                 !!! ---> PROBLEMS WITH ALLOCATED ARRAYS 
+                 !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
                  clearsky=.true.
                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
@@ -733,14 +736,8 @@
                       tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
                       zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
-                      fluxabs_sw1,fluxtop_dn,reffrad,tau_col1,             &
-                      cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
-
-                 clearsky=.false.
-                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
-                      albedo,emis,mu0,pplev,pplay,pt,                      &
-                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
-                      zdtlw2,zdtsw2,fluxsurf_lw2,fluxsurf_sw2,fluxtop_lw2, &
-                      fluxabs_sw2,fluxtop_dn,reffrad,tau_col2,             &
-                      cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
+                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
+                      reffrad,tau_col1,cloudfrac,totcloudfrac,             &
+		      clearsky,firstcall,lastcall)
+		 clearsky = .false.  !! just in case.     
 
                  ! Sum the fluxes and heating rates from cloudy/clear cases
@@ -749,28 +746,25 @@
                     ntf=1.-tf
                     
-                    fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw2(ig)
-                    fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw2(ig)
-                    fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw2(ig)
-                    fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw2(ig)
-                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col2(ig)
+                    fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
+                    fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
+                    fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
+                    fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
+                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
                     
                     do l=1,nlayer
-                       zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw2(ig,l)
-                       zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw2(ig,l)
+                       zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw(ig,l)
+                       zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw(ig,l)
                     enddo
+
+		    do nw=1,L_NSPECTV
+                       OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw)		       
+		    enddo
+		    do nw=1,L_NSPECTI
+                       OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw)		       
+		    enddo
+
                  enddo
 
-!             Otherwise standard callcorrk
-              else
-
-                 clearsky=.false.
-                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,         &
-                      albedo,emis,mu0,pplev,pplay,pt,                    &
-                      tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
-                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-                      fluxabs_sw,fluxtop_dn,reffrad,tau_col,             &
-                      cloudfrac,totcloudfrac,clearsky,firstcall,lastcall)
-
-              endif
+              endif !CLFvarying
 
 !             Radiative flux from the sky absorbed by the surface (W.m-2)
@@ -872,6 +866,8 @@
             print*,'In corrk SW atmospheric energy change =',dEtotSW,' W m-2'
             print*,'In corrk LW atmospheric energy change =',dEtotLW,' W m-2'
+            print*,'atmospheric energy change   (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
             print*,'In corrk SW surface energy change     =',dEtotsSW,' W m-2'
             print*,'In corrk LW surface energy change     =',dEtotsLW,' W m-2'
+            print*,'surface energy change       (SW+LW)   =',dEtotsLW+dEtotsSW,' W m-2'
          endif
 !-------------------------
@@ -884,5 +880,5 @@
 
       if (calldifv) then
-
+      
          do ig=1,ngrid
             zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
@@ -1743,4 +1739,5 @@
             OLR = OLR + area(ig)*fluxtop_lw(ig) 
             GND = GND + area(ig)*fluxgrd(ig) 
+	    if(.not.callsoil) GND=GND+ area(ig)*fluxrad(ig) 
             DYN = DYN + area(ig)*fluxdyn(ig) 
             
@@ -1765,8 +1762,9 @@
          if(enertest)then
             print*,'SW energy balance SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
+            print*,'LW energy balance LW++ + ***ASR*** = ',dEtotLW+dEtotsLW+ASR/totarea,' W m-2'
          endif
 
          if(meanOLR)then
-            if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
+            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
                ! to record global radiative balance
                open(92,file="rad_bal.out",form='formatted',access='append')
@@ -1857,5 +1855,5 @@
 
          if(meanOLR)then
-            if((ngridmx.gt.1) .or. (countG1D.eq.saveG1D))then
+            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
                ! to record global water balance
                open(98,file="h2o_bal.out",form='formatted',access='append')
@@ -1891,5 +1889,4 @@
 
 
-      if (ngrid.ne.1) then
          print*,''
          print*,'--> Ls =',zls*180./pi
@@ -1977,4 +1974,9 @@
 !                        "Solar radiative flux to space","W.m-2",2,         &
 !                        fluxtop_sw_tot)
+
+            call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
+            call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
+            call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
+
             call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
             call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
@@ -1984,4 +1986,15 @@
 
            if (tracer) then
+             do iq=1,nq
+                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
+                call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+                          'kg m^-2',2,qsurf(1,iq) )
+
+                call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+                          'kg m^-2',2,qcol(1,iq) )
+                call wstats(ngridmx,trim(noms(iq))//'_reff',                          & 
+                          trim(noms(iq))//'_reff',                                   & 
+                          'm',3,reffrad(1,1,iq))
+              end do
             if (water) then
                vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
@@ -2041,8 +2054,8 @@
 
 !     Output aerosols
-        call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
-        call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
-        call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
-        call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
+        if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
+        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
+        if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
+        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
 
 !     Output tracers
@@ -2093,55 +2106,10 @@
        endif
 
-!        ----------------------------------------------------------
-!        Output in netcdf file "diagsoil.nc" for subterranean
-!        variables (output every "ecritphy", as for writediagfi)
-!        ----------------------------------------------------------
-
-       ! Write soil temperature
-       !call writediagsoil(ngrid,"soiltemp","Soil temperature","K",3,tsoil)
-
-      else         ! if(ngrid.eq.1)
-
-!        ----------------------------------------------------------------------
-!        Output in grads file "g1d" (ONLY when using testphys1d)
-!        ----------------------------------------------------------------------
-
-         if(countG1D.eq.saveG1D)then
-         
-            call writeg1d(ngrid,1,tsurf,'tsurf','K')
-            call writeg1d(ngrid,1,ps,'ps','Pa')
-            call writeg1d(ngrid,nlayer,zt,'T','K')
-            call writeg1d(ngrid,nlayer,pq(1,1,1),'q','kg/kg')
-            call writeg1d(ngrid,nlayer,aerosol,'aerosol','SI')
-            call writeg1d(ngrid,nlayer,reffrad,'reffrad','SI')
-            call writeg1d(ngrid,nlayer,zdtlw,'dtlw','SI')
-            call writeg1d(ngrid,nlayer,zdtsw,'dtsw','SI')
-            call writeg1d(ngrid,nlayer,zdtdyn,'dtdyn','SI')
-
-!           radiation balance
-            call writeg1d(ngrid,1,ISR/totarea,'ISR','W m-2')
-            call writeg1d(ngrid,1,ASR/totarea,'ASR','W m-2')
-            call writeg1d(ngrid,1,OLR/totarea,'OLR','W m-2')  
-
-            if(tracer) then
-               do iq=1,nq
-                  call writeg1d(ngrid,1,qsurf(1,iq),trim(noms(iq))//'_s','kg m^-2') 
-                  call writeg1d(ngrid,1,qcol(1,iq),trim(noms(iq))//'_c','kg m^-2') 
-                  call writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 
-               end do
-               
-               if(waterrain)then
-                  CALL writeg1d(ngrid,1,zdqsrain,'rainfall','kg m-2 s-1')
-                  CALL writeg1d(ngrid,1,zdqssnow,'snowfall','kg m-2 s-1')
-               endif
-               
-            end if
-
-            countG1D=1
-         else
-            countG1D=countG1D+1
-         endif ! if time to save
-
-      endif        ! if(ngrid.ne.1)
+!      output spectrum 
+      if(specOLR.and.corrk)then      
+                 call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
+                 call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
+      endif
+
 
       icount=icount+1
Index: /trunk/LMDZ.GENERIC/libf/phystd/rain.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 526)
@@ -62,5 +62,5 @@
 
       logical simple                    ! Use very simple Emanuel scheme?
-      parameter(simple=.true.)
+      parameter(simple=.false.)
 
       logical evap_prec                 ! Does the rain evaporate?
Index: /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 526)
@@ -41,5 +41,4 @@
 #include "comvert.h"
 #include "netcdf.inc"
-#include "comg1d.h"
 #include "watercap.h" 
 #include "fisice.h"
@@ -117,5 +116,5 @@
 c=======================================================================
 
-      saveprofile=.true.
+      saveprofile=.false.
 
 c ------------------------------------------------------
@@ -274,4 +273,10 @@
       write(*,*)" time = ",time
       time=time/24.E+0 ! convert time (hours) to fraction of sol
+      
+       ecritphy=day_step ! default value for ecritphy
+       PRINT *,'Nunber of steps between writediagfi ?'
+       call getin("ecritphy",ecritphy)
+       write(*,*) " ecritphy = ",ecritphy
+
 
 c  Discretisation (Definition de la grille et des pas de temps)
@@ -294,4 +299,11 @@
       ndt=ndt*day_step     
       dtphys=daysec/day_step  
+
+
+c output spectrum?
+      write(*,*)"Output spectral OLR?"
+      specOLR=.false.
+      call getin("specOLR",specOLR)
+      write(*,*)" specOLR = ",specOLR
 
 
@@ -304,6 +316,6 @@
       write(*,*) " psurf = ",psurf
 c Pression de reference
-      pa=20.
-      preff=610. ! these values are not needed in 1D  
+      pa=psurf/30.
+      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
 
 c Gravity of planet
@@ -577,16 +589,4 @@
       enddo
 
-c    Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
-c    ----------------------------------------------------------------
-c    (effectuee avec "writeg1d" a partir de "physiq.F"
-c    ou tout autre subroutine physique, y compris celle ci).
-
-        g1d_nlayer=nlayer
-        g1d_nomfich='g1d.dat'
-        g1d_unitfich=40
-        g1d_nomctl='g1d.ctl'
-        g1d_unitctl=41
-        g1d_premier=.true.
-        g2d_premier=.true.
 
 c  Ecriture de "startfi"
@@ -613,13 +613,13 @@
          lastcall=.true.
          call stellarlong(day*1.0,zls)
-         write(103,*) 'Ls=',zls*180./pi
-         write(103,*) 'Lat=', lati(1)*180./pi
-         write(103,*) 'Tau=', tauvis/700*psurf
-         write(103,*) 'RunEnd - Atmos. Temp. File'
-         write(103,*) 'RunEnd - Atmos. Temp. File'
-         write(104,*) 'Ls=',zls*180./pi
-         write(104,*) 'Lat=', lati(1)
-         write(104,*) 'Tau=', tauvis/700*psurf
-         write(104,*) 'RunEnd - Atmos. Temp. File'
+!         write(103,*) 'Ls=',zls*180./pi
+!         write(103,*) 'Lat=', lati(1)*180./pi
+!         write(103,*) 'Tau=', tauvis/700*psurf
+!         write(103,*) 'RunEnd - Atmos. Temp. File'
+!         write(103,*) 'RunEnd - Atmos. Temp. File'
+!         write(104,*) 'Ls=',zls*180./pi
+!         write(104,*) 'Lat=', lati(1)
+!         write(104,*) 'Tau=', tauvis/700*psurf
+!         write(104,*) 'RunEnd - Atmos. Temp. File'
         ENDIF
 
@@ -748,8 +748,4 @@
 c    ========================================================
 
-c    fermeture pour conclure les sorties format grads dans "g1d.dat"
-c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
-c    ou tout autre subroutine physique, y compris celle ci).
-
 ! save temperature profile
       if(saveprofile)then
@@ -761,9 +757,4 @@
       endif
 
-      logplevs(1)=log(plev(1)/100000.0)
-      do ilayer=2,nlayer
-         logplevs(ilayer)=log(play(ilayer)/100000.0)
-      enddo
-      CALL endg1d(1,nlayer,logplevs,ndt) ! now simply log(p) with p in bars
 
 c    ========================================================
Index: /trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F	(revision 526)
@@ -64,5 +64,5 @@
 C     TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
 
-      TTOP  = TLEV(2)
+      TTOP  = TLEV(2)  ! JL12 why not (1) ???
       TSURF = TLEV(L_LEVELS)
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F	(revision 526)
@@ -196,5 +196,5 @@
         ! are undefined; so set them to 1
         iphysiq=1 
-        irythme=1
+        !JL11 irythme=1
         ! NB: 
       endif
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 526)
@@ -112,5 +112,6 @@
 ! The following test is here to enforce that writediagfi is not used with the
 ! 1D version of the GCM
-      if (ngrid.eq.1) return
+!not anymore (JL12)
+      if (ngrid.eq.-1) return
       
 c     nom=trim((nom))
@@ -238,5 +239,5 @@
               ierr= NF_INQ_DIMID(nid,"longitude",id(1))
               ierr= NF_INQ_DIMID(nid,"latitude",id(2))
-              ierr= NF_INQ_DIMID(nid,"bandsIR",id(3))
+              ierr= NF_INQ_DIMID(nid,"IR Wavenumber",id(3))
               ierr= NF_INQ_DIMID(nid,"Time",id(4))
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 525)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 526)
@@ -112,5 +112,6 @@
 ! The following test is here to enforce that writediagfi is not used with the
 ! 1D version of the GCM
-      if (ngrid.eq.1) return
+!not anymore (JL12)
+      if (ngrid.eq.-1) return
       
 c     nom=trim((nom))
@@ -238,5 +239,5 @@
               ierr= NF_INQ_DIMID(nid,"longitude",id(1))
               ierr= NF_INQ_DIMID(nid,"latitude",id(2))
-              ierr= NF_INQ_DIMID(nid,"bandsVI",id(3))
+              ierr= NF_INQ_DIMID(nid,"VI Wavenumber",id(3))
               ierr= NF_INQ_DIMID(nid,"Time",id(4))
 
