Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 961)
@@ -113,5 +113,5 @@
 
       REAL*8 tauaero(L_LEVELS+1,naerkind)
-      REAL*8 nfluxtopv,nfluxtopi,nfluxtop
+      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
       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)
@@ -125,5 +125,4 @@
       INTEGER icount
 
-      real fluxtoplanet
       real szangle
       logical global1d
@@ -291,8 +290,6 @@
 
 !     how much light we get
-      fluxtoplanet=0
       do nw=1,L_NSPECTV
          stel(nw)=stellarf(nw)/(dist_star**2)
-         fluxtoplanet=fluxtoplanet + stel(nw)
       end do
 
@@ -657,10 +654,7 @@
          if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
 
-            fluxtoplanet=0.
-
             if((ngrid.eq.1).and.(global1d))then
                do nw=1,L_NSPECTV
                   stel_fract(nw)= stel(nw) * 0.25 / acosz
-                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
                                 ! globally averaged = divide by 4
                                 ! but we correct for solar zenith angle
@@ -669,5 +663,4 @@
                do nw=1,L_NSPECTV
                   stel_fract(nw)= stel(nw) * fract(ig)
-                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
                end do
             endif
@@ -679,5 +672,5 @@
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
                  acosz,stel_fract,gweight,                         &
-                 nfluxtopv,nfluxoutv_nu,nfluxgndv_nu,              &
+                 nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,              &
                  fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
 
@@ -709,5 +702,5 @@
 
 !     Flux incident at the top of the atmosphere
-         fluxtop_dn(ig)=fluxdnv(1)
+         fluxtop_dn(ig)=fluxtopvdn 
 
          fluxtop_lw(ig)  = real(nfluxtopi)
Index: trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F	(revision 961)
@@ -156,5 +156,5 @@
         EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL 
         EM    = 1.0D0/EP
-        TERM  = UBARI/(1.-W0(L)*COSBAR(L))
+        TERM  = UBARI/(1.D0-W0(L)*COSBAR(L))
 
 C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 961)
@@ -549,6 +549,6 @@
                write(*,*) "cpp=",cpp
            ENDIF
-         else
-           mugaz=8.314*1000./pr
+!         else
+!           mugaz=8.314*1000./pr
          endif
          call su_gases
Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 961)
@@ -120,5 +120,5 @@
      ! if we have continuum opacities, we need dz
      if(kastprof)then
-        dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K))
+        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
         U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
      else
@@ -149,5 +149,5 @@
 
 ! continuum absorption
-        DCONT = 0.0 
+        DCONT = 0.0d0 
 
         if(continuum.and.(.not.graybody))then
@@ -163,5 +163,5 @@
               endif
 
-              dtemp=0.0
+              dtemp=0.0d0
               if(igas.eq.igas_N2)then
 
@@ -180,5 +180,5 @@
                  do jgas=1,ngasmx
                     p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
-                    dtempc  = 0.0
+                    dtempc  = 0.0d0
                     if(jgas.eq.igas_N2)then 
                        interm = indi(nw,igas,jgas)
@@ -286,5 +286,5 @@
   end do
 
-
+  DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0
   !=======================================================================
   !     Now the full treatment for the layers, where besides the opacity
@@ -315,5 +315,5 @@
 
         atemp = 0.
-        if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
+        if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
            do iaer=1,naerkind
               atemp = atemp +                                     &
@@ -324,8 +324,8 @@
         else
            WBARI(L,nw,ng) = 0.0D0
-           DTAUI(L,NW,NG) = 1.0E-9
+           DTAUI(L,NW,NG) = 1.0D-9
         endif
 
-        if(btemp(L,nw) .GT. 0.0) then
+        if(btemp(L,nw) .GT. 0.0d0) then
            cosbi(L,NW,NG) = atemp/btemp(L,nw)
         else
@@ -344,5 +344,5 @@
               DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
 
-              if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
+              if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
 
                  WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
@@ -350,5 +350,5 @@
               else
                  WBARI(L,nw,ng) = 0.0D0
-                 DTAUI(L,NW,NG) = 1.0E-9
+                 DTAUI(L,NW,NG) = 1.0D-9
               endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 961)
@@ -1433,4 +1433,5 @@
          print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
          print*, ISR,ASR,OLR,GND,DYN
+
          if(enertest)then
             print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
Index: trunk/LMDZ.GENERIC/libf/phystd/setspi.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/setspi.F90	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/setspi.F90	(revision 961)
@@ -127,10 +127,10 @@
 
       do M=1,L_NSPECTI
-         WNOI(M)  = 0.5*(BWNI(M+1)+BWNI(M))
+         WNOI(M)  = 0.5D0*(BWNI(M+1)+BWNI(M))
          DWNI(M)  = BWNI(M+1)-BWNI(M)
-         WAVEI(M) = 1.0E+4/WNOI(M)
-         BLAMI(M) = 0.01/BWNI(M)         
-      end do
-      BLAMI(M) = 0.01/BWNI(M)
+         WAVEI(M) = 1.0D+4/WNOI(M)
+         BLAMI(M) = 0.01D0/BWNI(M)         
+      end do
+      BLAMI(M) = 0.01D0/BWNI(M)
 !     note M=L_NSPECTI+1 after loop due to Fortran bizarreness
 
@@ -148,6 +148,6 @@
          a = 1.0D-2/BWNI(NW+1)
          b = 1.0D-2/BWNI(NW)
-         bpa = (b+a)/2.0
-         bma = (b-a)/2.0
+         bpa = (b+a)/2.0D0
+         bma = (b-a)/2.0D0
          do nt=NTstar,NTstop
             T   = dble(NT)/NTfac
@@ -167,5 +167,5 @@
          print*,'setspi: Force F=sigma*eps*T^4 for all values of T!'
          do nt=NTstar,NTstop
-            plancksum=0.0
+            plancksum=0.0D0
             T=dble(NT)/NTfac
        
@@ -185,5 +185,5 @@
       if(planckcheck)then
          ! check energy conservation at lower temperature boundary
-         plancksum=0.0
+         plancksum=0.0D0
          nt=NTstar
          do NW=1,L_NSPECTI
@@ -195,5 +195,5 @@
          
          ! check energy conservation at upper temperature boundary
-         plancksum=0.0
+         plancksum=0.0D0
          nt=NTstop
          do NW=1,L_NSPECTI
Index: trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F	(revision 961)
@@ -70,6 +70,4 @@
       NTS   = int(TSURF*NTfac)-NTstar+1
       NTT   = int(TTOP *NTfac)-NTstar+1
-!      NTS   = TSURF*10.0D0-499
-!      NTT   = TTOP *10.0D0-499
 
 !JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4
@@ -110,5 +108,6 @@
 C         OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
  
-          TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+!          TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+          TAUTOP = TAUCUMI(2,NW,NG)
           BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
  
@@ -161,5 +160,6 @@
        NG     = L_NGAUSS
 
-       TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+!       TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
+       TAUTOP = TAUCUMI(2,NW,NG)
        BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
 
Index: trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F	(revision 961)
@@ -1,5 +1,5 @@
       SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
-     *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUXOUTV_nu,
-     *                  NFLUXGNDV_nu,
+     *                  UBAR0,STEL,GWEIGHT,NFLUXTOPV,FLUXTOPVDN,
+     *                  NFLUXOUTV_nu,NFLUXGNDV_nu,
      *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
 
@@ -18,5 +18,5 @@
       real*8 STEL(L_NSPECTV)
       real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
-      real*8 NFLUXTOPV, FLUXUP, FLUXDN
+      real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN
       real*8 NFLUXOUTV_nu(L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
@@ -37,4 +37,5 @@
 
       NFLUXTOPV = 0.0
+      FLUXTOPVDN = 0.0
 
       DO NW=1,L_NSPECTV
@@ -100,4 +101,6 @@
           NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*GWEIGHT(NG)*
      *                          (1.0-FZEROV(NW))
+          FLUXTOPVDN = FLUXTOPVDN+FLUXDN*GWEIGHT(NG)*
+     *                          (1.0-FZEROV(NW))
           DO L=1,L_NLAYRAD
             FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*
@@ -156,4 +159,5 @@
 
         NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*FZERO
+        FLUXTOPVDN = FLUXTOPVDN+FLUXDN*FZERO
         DO L=1,L_NLAYRAD
           FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*FZERO
Index: trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 959)
+++ trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 961)
@@ -372,5 +372,5 @@
 
          do nw=1,L_NSPECTI
-            fzeroi(nw) = 0.
+            fzeroi(nw) = 0.d0
 !            do nt=1,L_NTREF
 !               do np=1,L_NPREF
@@ -378,5 +378,5 @@
 !                     do ng = 1,L_NGAUSS
 !                        if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzeroi(nw)=fzeroi(nw)+1.
+!                           fzeroi(nw)=fzeroi(nw)+1.d0
 !                        endif
 !                     end do
@@ -388,5 +388,5 @@
 
          do nw=1,L_NSPECTV
-            fzerov(nw) = 0.
+            fzerov(nw) = 0.d0
 !            do nt=1,L_NTREF
 !               do np=1,L_NPREF
@@ -394,5 +394,5 @@
 !                     do ng = 1,L_NGAUSS
 !                        if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzerov(nw)=fzerov(nw)+1.
+!                           fzerov(nw)=fzerov(nw)+1.d0
 !                        endif
 !                     end do
