Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 486)
+++ 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
