Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 955)
+++ trunk/LMDZ.GENERIC/README	(revision 959)
@@ -937,2 +937,5 @@
 - moved the 'makbands' script from "grid" to (more appropriate) "phystd/bands"
   subdirectory, and consequently adapted the makegcm_* scripts
+
+== 15/05/2013 == JL
+- correction in radiative scheme to enforce double precision
Index: trunk/LMDZ.GENERIC/libf/phystd/blackl.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/blackl.F	(revision 955)
+++ trunk/LMDZ.GENERIC/libf/phystd/blackl.F	(revision 959)
@@ -4,5 +4,5 @@
 
       ! physical constants
-      sigma=5.6693d-08
+      sigma=5.67032D-8
       pi=datan(1.d0)*4.d0
       c0=2.9979d+08
Index: trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F	(revision 955)
+++ trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F	(revision 959)
@@ -75,14 +75,12 @@
 !           endif
 ! Prevent this with an if statement:
-        if (W0(L).eq.1.) then
-           W0(L) = 0.99999
+        if (W0(L).eq.1.D0) then
+           W0(L) = 0.99999D0
         endif
 !----------------------------------------------------
 
-        ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) )
-        LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI
-
-        !NT2   = TLEV(2*L+2)*10.0D0-499
-        !NT    = TLEV(2*L)*10.0D0-499
+        ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
+        LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
+
         NT    = int(TLEV(2*L)*NTfac)   - NTstar+1
         NT2   = int(TLEV(2*L+2)*NTfac) - NTstar+1
@@ -97,12 +95,10 @@
 
       if (W0(L).eq.1.) then
-          W0(L) = 0.99999
+          W0(L) = 0.99999D0
       end if
 
-      ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) )
-      LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI
-
-      !NT    = TLEV(2*L+1)*10.0D0-499
-      !NT2   = TLEV(2*L)*10.0D0-499
+      ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
+      LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
+
       NT    = int(TLEV(2*L+1)*NTfac) - NTstar+1
       NT2   = int(TLEV(2*L)*NTfac)   - NTstar+1
@@ -111,6 +107,6 @@
 
       DO L=1,L_NLAYRAD
-        GAMA(L) = (1.0-ALPHA(L))/(1.0+ALPHA(L))
-        TERM    = UBARI/(1.0-W0(L)*COSBAR(L))
+        GAMA(L) = (1.0D0-ALPHA(L))/(1.0D0+ALPHA(L))
+        TERM    = UBARI/(1.0D0-W0(L)*COSBAR(L))
 
 C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
@@ -137,5 +133,5 @@
 
         EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX))
-        EM    = 1.0/EP
+        EM    = 1.0D0/EP
         E1(L) = EP+GAMA(L)*EM
         E2(L) = EP-GAMA(L)*EM
@@ -159,5 +155,5 @@
         DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L)
         EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL 
-        EM    = 1.0/EP
+        EM    = 1.0D0/EP
         TERM  = UBARI/(1.-W0(L)*COSBAR(L))
 
@@ -181,6 +177,6 @@
 
       EP   = EXP(MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL 
-      EM   = 1.0/EP
-      TERM = UBARI/(1.-W0(L)*COSBAR(L))
+      EM   = 1.0D0/EP
+      TERM = UBARI/(1.D0-W0(L)*COSBAR(L))
 
 C     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
@@ -199,7 +195,7 @@
 C     FLUX AT THE PTOP LEVEL
 
-      EP   = 1.0
-      EM   = 1.0
-      TERM = UBARI/(1.0-W0(1)*COSBAR(1))
+      EP   = 1.0D0
+      EM   = 1.0D0
+      TERM = UBARI/(1.0D0-W0(1)*COSBAR(1))
 
 C     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 955)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 959)
@@ -326,5 +326,5 @@
       real mass(ngrid,nlayermx),massarea(ngrid,nlayermx)
       real dEtot, dEtots, AtmToSurf_TurbFlux
-      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
+      real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
       real dEzRadsw(ngrid,nlayermx),dEzRadlw(ngrid,nlayermx),dEzdiff(ngrid,nlayermx)
       real dEdiffs(ngrid),dEdiff(ngrid)
@@ -816,5 +816,5 @@
             dEtotSW  = cpp*SUM(massarea(:,:)*zdtsw(:,:))/totarea
             dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
-            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
+            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea !JL13 carefull, albedo can have changed since the last time we called corrk
             dEtotsLW = SUM((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:))/totarea
 	    dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
@@ -1189,5 +1189,5 @@
            ! find qtot
            if(watertest)then
-              iq=3
+              iq=igcm_h2o_ice
 	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
 	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
@@ -1203,5 +1203,5 @@
            ! find qtot
            if(watertest)then
-              iq=3
+              iq=igcm_h2o_ice
 	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
 	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
@@ -1397,5 +1397,10 @@
            print*,Ts1,Ts2,Ts3,TsS
          endif
-      endif
+      else
+         if (cursor == 1) then
+           print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
+           print*,Ts1,Ts2,Ts3
+         endif         
+      end if
 
 !     ---------------------------------------------------------
@@ -1426,7 +1431,7 @@
          endif
 
+         print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
+         print*, ISR,ASR,OLR,GND,DYN
          if(enertest)then
-            print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
-            print*, ISR,ASR,OLR,GND,DYN
             print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
             print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
@@ -1765,6 +1770,6 @@
           do ig=1,ngrid
              if(tau_col(ig).gt.1.e3)then
-                print*,'WARNING: tau_col=',tau_col(ig)
-                print*,'at ig=',ig,'in PHYSIQ'
+!                print*,'WARNING: tau_col=',tau_col(ig)
+!                print*,'at ig=',ig,'in PHYSIQ'
              endif
           end do
Index: trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 955)
+++ trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 959)
@@ -128,5 +128,5 @@
       real*8, parameter :: SCALEP = 1.00D+2
 
-      real*8, parameter :: sigma = 5.67032e-8
+      real*8, parameter :: sigma = 5.67032D-8
 
       real*8 Cmk
Index: trunk/LMDZ.GENERIC/libf/phystd/setspi.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/setspi.F90	(revision 955)
+++ trunk/LMDZ.GENERIC/libf/phystd/setspi.F90	(revision 959)
@@ -68,5 +68,5 @@
       mm=0
 
-      forceEC=.false.
+      forceEC=.true.
       planckcheck=.true.
 
Index: trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F	(revision 955)
+++ trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F	(revision 959)
@@ -43,13 +43,13 @@
 C     ZERO THE NET FLUXES
     
-      NFLUXTOPI = 0.0
+      NFLUXTOPI = 0.0D0
 
       DO NW=1,L_NSPECTI
-        NFLUXTOPI_nu(NW) = 0.0
+        NFLUXTOPI_nu(NW) = 0.0D0
         DO L=1,L_NLAYRAD
-           FLUXUPI_nu(L,NW) = 0.0
-
-              fup_tmp(nw)=0.0
-              fdn_tmp(nw)=0.0
+           FLUXUPI_nu(L,NW) = 0.0D0
+
+              fup_tmp(nw)=0.0D0
+              fdn_tmp(nw)=0.0D0
 
         END DO
@@ -57,7 +57,7 @@
 
       DO L=1,L_NLAYRAD
-        FMNETI(L)  = 0.0
-        FLUXUPI(L) = 0.0
-        FLUXDNI(L) = 0.0
+        FMNETI(L)  = 0.0D0
+        FLUXUPI(L) = 0.0D0
+        FLUXDNI(L) = 0.0D0
       END DO
  
@@ -102,5 +102,5 @@
          
           if(TAUGSURF(NW,NG).lt. TLIMIT) then
-            fzero = fzero + (1.0-FZEROI(NW))*GWEIGHT(NG)
+            fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG)
             goto 30
           end if
@@ -111,5 +111,5 @@
  
           TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
-          BTOP   = (1.0-EXP(-TAUTOP/UBARI))*PLTOP
+          BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
  
 C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
@@ -127,9 +127,9 @@
 
           NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG)*
-     *                           (1.0-FZEROI(NW))
+     *                           (1.0D0-FZEROI(NW))
 
 c         and same thing by spectral band... (RDW)
           NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)
-     *      +FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))
+     *      +FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
 
 
@@ -139,13 +139,13 @@
 
             FMNETI(L)  = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*
-     *                              GWEIGHT(NG)*(1.0-FZEROI(NW))
+     *                              GWEIGHT(NG)*(1.0D0-FZEROI(NW))
             FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*
-     *                                (1.0-FZEROI(NW))
+     *                                (1.0D0-FZEROI(NW))
             FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*GWEIGHT(NG)*
-     *                                (1.0-FZEROI(NW))
+     *                                (1.0D0-FZEROI(NW))
 
 c         and same thing by spectral band... (RW)
             FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + 
-     *                FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW))
+     *                FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0D0-FZEROI(NW))
 
           END DO
@@ -162,5 +162,5 @@
 
        TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2))
-       BTOP   = (1.0-EXP(-TAUTOP/UBARI))*PLTOP
+       BTOP   = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP
 
 C      WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
