Index: LMDZ6/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 3748)
+++ LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 3756)
@@ -142,4 +142,5 @@
         <field id="SWdnSFCclr"    long_name="SWdn clear sky at surface"    unit="W/m2" />
         <field id="SWdnSFCcleanclr"    long_name="SWdn clean (no aerosol) clear sky at surface"    unit="W/m2" />
+        <field id="fdiffSWdnSFC"    long_name="Fraction of diffuse SWdn at surface"    unit="-" />
         <field id="LWupSFC"    long_name="Upwd. IR rad. at surface"    unit="W/m2" />
         <field id="LWupSFCclr"    long_name="CS Upwd. IR rad. at surface"    unit="W/m2" />
Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3756)
@@ -164,5 +164,5 @@
        rlon,      rlat,      rugoro,   rmu0,          &
        zsig,      lwdown_m,  pphi,     cldt,          &
-       rain_f,    snow_f,    solsw_m,  sollw_m,       &
+       rain_f,    snow_f,    solsw_m,  solswfdiff_m, sollw_m,       &
        gustiness,                                     &
        t,         q,         u,        v,             &
@@ -318,4 +318,5 @@
     REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
     REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
+    REAL, DIMENSION(klon),        INTENT(IN)        :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface
     REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! net longwave radiation at mean surface
     REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t       ! temperature (K)
@@ -453,10 +454,10 @@
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxv    ! v wind tension, mean for each grid point
-    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)     :: z0m,z0h      ! rugosity length (m)
-    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: agesno   ! age of snow at surface
+    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: z0m,z0h      ! rugosity length (m)
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: agesno   ! age of snow at surface
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: solsw      ! net shortwave radiation at surface 
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw      ! net longwave radiation at surface
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts       ! change in temperature at surface
-    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: evap     ! evaporation at surface
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap     ! evaporation at surface
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat    ! latent flux
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m        ! temperature at 2 meter height
@@ -965,5 +966,5 @@
     yqsurf = 0.0  ; yalb = 0.0 ; yalb_vis = 0.0
 !albedo SB <<<
-    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0    
+    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0
     ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yu1 = 0.0    
     yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
@@ -1183,4 +1184,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+!--OB this line is not satisfactory because alb is the direct albedo not total albedo
           solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
        ENDDO
@@ -1205,4 +1207,8 @@
 !>al1
 
+!--OB add diffuse fraction of SW down
+   DO n=1,nbcf_out
+       IF (cfname_out(n) == "swdownfdiff" ) fields_out(:,n) = solswfdiff_m(:)
+   ENDDO
 ! >> PC
    IF (carbon_cycle_cpl .AND. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
Index: LMDZ6/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 3756)
@@ -15,5 +15,5 @@
        ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, radpas, radsol, rain_fall, ratqs, &
        rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
-       solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
+       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
        wake_s, wake_dens, zgam, zmax0, zmea, zpic, zsig, &
@@ -307,11 +307,12 @@
 
   found=phyetat0_get(1,solsw,"solsw","net SW radiation surf",0.)
+  found=phyetat0_get(1,solswfdiff,"solswfdiff","fraction of SW radiation surf that is diffuse",1.)
   found=phyetat0_get(1,sollw,"sollw","net LW radiation surf",0.)
   found=phyetat0_get(1,sollwdown,"sollwdown","down LW radiation surf",0.)
   IF (.NOT. found) THEN
-     sollwdown = 0. ;  zts=0.
-     do nsrf=1,nbsrf
+     sollwdown(:) = 0. ;  zts(:)=0.
+     DO nsrf=1,nbsrf
         zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
-     enddo
+     ENDDO
      sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
   ENDIF
Index: LMDZ6/trunk/libf/phylmd/phyredem.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 3756)
@@ -27,5 +27,5 @@
                                 ale_wake, ale_bl_stat,                       &
                                 du_gwd_rando, du_gwd_front, u10m, v10m,      &
-                                treedrg
+                                treedrg, solswfdiff
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
   USE iostart, ONLY: open_restartphy, close_restartphy, enddef_restartphy, put_field, put_var
@@ -184,4 +184,6 @@
 
     CALL put_field(pass,"solsw", "Rayonnement solaire a la surface", solsw)
+
+    CALL put_field(pass,"solswfdiff", "Fraction du rayonnement solaire a la surface qui est diffus", solswfdiff)
 
     CALL put_field(pass,"sollw", "Rayonnement IF a la surface", sollw)
Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3756)
@@ -461,5 +461,7 @@
   TYPE(ctrl_out), SAVE :: o_SWupSFCcleanclr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'SWupSFCcleanclr', 'SWup clear sky clean (no aerosol) at surface', 'W/m2', (/ ('', i=1, 10) /))
-  TYPE(ctrl_out), SAVE :: o_SWdnSFC = ctrl_out((/ 1, 1, 10, 10, 5, 10, 11, 11, 11, 11/), &
+  TYPE(ctrl_out), SAVE :: o_fdiffSWdnSFC = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'fdiffSWdnSFC', 'Fraction of diffuse SWdn at surface', 'W/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_SWdnSFC = ctrl_out((/ 1, 1, 1, 10, 5, 10, 11, 11, 11, 11/), &
     'SWdnSFC', 'SWdn at surface', 'W/m2', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_SWdnSFCclr = ctrl_out((/ 1, 4, 10, 10, 5, 10, 11, 11, 11, 11/), &
Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3756)
@@ -45,5 +45,5 @@
          o_snow, o_msnow, o_fsnow, o_evap, o_ep,o_epmax_diag, & ! epmax_cape
          o_tops, o_tops0, o_topl, o_topl0, &
-         o_SWupTOA, o_SWupTOAclr, o_SWupTOAcleanclr, o_SWdnTOA, &
+         o_SWupTOA, o_SWupTOAclr, o_SWupTOAcleanclr, o_SWdnTOA, o_fdiffSWdnSFC, &
          o_SWdnTOAclr, o_nettop, o_SWup200, &
          o_SWup200clr, o_SWdn200, o_SWdn200clr, &
@@ -223,5 +223,5 @@
          qsol, z0m, z0h, fevap, agesno, &
          nday_rain, rain_con, snow_con, &
-         topsw, toplw, toplw0, swup, swdn, &
+         topsw, toplw, toplw0, swup, swdn, solswfdiff, &
          topsw0, swupc0, swdnc0, swup0, swdn0, SWup200, SWup200clr, &
          SWdn200, SWdn200clr, LWup200, LWup200clr, &
@@ -938,4 +938,6 @@
        ENDIF
        CALL histwrite_phy(o_SWdnSFCcleanclr, zx_tmp_fi2d)
+
+       CALL histwrite_phy(o_fdiffSWdnSFC, solswfdiff)
 
        IF (vars_defined) THEN
Index: LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3756)
@@ -47,5 +47,5 @@
 !albedo SB >>>
       REAL, ALLOCATABLE, SAVE :: falb_dif(:,:,:), falb_dir(:,:,:)
-      real, allocatable, save :: chl_con(:)
+      REAL, ALLOCATABLE, SAVE :: chl_con(:)
 !$OMP THREADPRIVATE(falb_dir,falb_dif,chl_con)
 !albedo SB <<<
@@ -54,6 +54,6 @@
       REAL, ALLOCATABLE, SAVE :: rain_fall(:), snow_fall(:)
 !$OMP THREADPRIVATE( rain_fall, snow_fall)
-      REAL, ALLOCATABLE, SAVE :: solsw(:), sollw(:)
-!$OMP THREADPRIVATE(solsw, sollw)
+      REAL, ALLOCATABLE, SAVE :: solsw(:), solswfdiff(:), sollw(:)
+!$OMP THREADPRIVATE(solsw, solswfdiff, sollw)
       REAL, ALLOCATABLE, SAVE :: radsol(:)
 !$OMP THREADPRIVATE(radsol)
@@ -358,6 +358,6 @@
       REAL,ALLOCATABLE,SAVE :: albplap(:)
 !$OMP THREADPRIVATE(albplap)
-      REAL,ALLOCATABLE,SAVE :: solswp(:), sollwp(:)
-!$OMP THREADPRIVATE(solswp, sollwp)
+      REAL,ALLOCATABLE,SAVE :: solswp(:), solswfdiffp(:), sollwp(:)
+!$OMP THREADPRIVATE(solswp, solswfdiffp, sollwp)
       REAL,ALLOCATABLE,SAVE :: sollwdownp(:)
 !$OMP THREADPRIVATE(sollwdownp)
@@ -466,5 +466,5 @@
       ALLOCATE(rain_fall(klon))
       ALLOCATE(snow_fall(klon))
-      ALLOCATE(solsw(klon), sollw(klon))
+      ALLOCATE(solsw(klon), solswfdiff(klon), sollw(klon))
       sollw=0.0
       ALLOCATE(radsol(klon))
@@ -597,5 +597,5 @@
       ALLOCATE(radsolp(klon), topswp(klon), toplwp(klon))
       ALLOCATE(albplap(klon))
-      ALLOCATE(solswp(klon), sollwp(klon))
+      ALLOCATE(solswp(klon), solswfdiffp(klon), sollwp(klon))
       ALLOCATE(gustiness(klon))
       ALLOCATE(sollwdownp(klon))
@@ -628,5 +628,5 @@
 !!! fin nrlmd le 10/04/2012
       IF (ok_gwd_rando) THEN
-        allocate(du_gwd_rando(klon, klev))
+        ALLOCATE(du_gwd_rando(klon, klev))
         du_gwd_rando(:,:)=0.
       ENDIF
@@ -644,127 +644,127 @@
 include "clesphys.h"
 
-      deallocate(pctsrf, ftsol, falb1, falb2)
-      deallocate(qsol,fevap,z0m,z0h,agesno)
+      DEALLOCATE(pctsrf, ftsol, falb1, falb2)
+      DEALLOCATE(qsol,fevap,z0m,z0h,agesno)
 !FC
-      deallocate(treedrg)
-      deallocate(rain_fall, snow_fall, solsw, sollw, radsol, swradcorr)
-      deallocate(zmea, zstd, zsig, zgam)
-      deallocate(zthe, zpic, zval)
-      deallocate(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
-      deallocate(qs_ancien, ql_ancien)
-      deallocate(prw_ancien, prlw_ancien, prsw_ancien)
-      deallocate(qtc_cv,sigt_cv)
-      deallocate(u_ancien, v_ancien)
-      deallocate(tr_ancien)                           !RomP
-      deallocate(ratqs, pbl_tke,coefh,coefm)
+      DEALLOCATE(treedrg)
+      DEALLOCATE(rain_fall, snow_fall, solsw, solswfdiff, sollw, radsol, swradcorr)
+      DEALLOCATE(zmea, zstd, zsig, zgam)
+      DEALLOCATE(zthe, zpic, zval)
+      DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
+      DEALLOCATE(qs_ancien, ql_ancien)
+      DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien)
+      DEALLOCATE(qtc_cv,sigt_cv)
+      DEALLOCATE(u_ancien, v_ancien)
+      DEALLOCATE(tr_ancien)                           !RomP
+      DEALLOCATE(ratqs, pbl_tke,coefh,coefm)
 !nrlmd<
-      deallocate(delta_tsurf)
+      DEALLOCATE(delta_tsurf)
 !>nrlmd
-      deallocate(zmax0, f0)
-      deallocate(sig1, w01)
-      deallocate(entr_therm, fm_therm)
-      deallocate(detr_therm)
-      deallocate(clwcon0th, rnebcon0th)
+      DEALLOCATE(zmax0, f0)
+      DEALLOCATE(sig1, w01)
+      DEALLOCATE(entr_therm, fm_therm)
+      DEALLOCATE(detr_therm)
+      DEALLOCATE(clwcon0th, rnebcon0th)
 ! radiation outputs
-      deallocate(swdnc0, swdn0, swdn)
-      deallocate(swupc0, swup0, swup)
-      deallocate(lwdnc0, lwdn0, lwdn)
-      deallocate(lwupc0, lwup0, lwup)
-      deallocate(SWdn200clr, SWdn200)
-      deallocate(SWup200clr, SWup200)
-      deallocate(LWdn200clr, LWdn200)
-      deallocate(LWup200clr, LWup200)
-      deallocate(LWdnTOA, LWdnTOAclr)
+      DEALLOCATE(swdnc0, swdn0, swdn)
+      DEALLOCATE(swupc0, swup0, swup)
+      DEALLOCATE(lwdnc0, lwdn0, lwdn)
+      DEALLOCATE(lwupc0, lwup0, lwup)
+      DEALLOCATE(SWdn200clr, SWdn200)
+      DEALLOCATE(SWup200clr, SWup200)
+      DEALLOCATE(LWdn200clr, LWdn200)
+      DEALLOCATE(LWup200clr, LWup200)
+      DEALLOCATE(LWdnTOA, LWdnTOAclr)
 ! pressure level
-      deallocate(tsumSTD)
-      deallocate(usumSTD, vsumSTD)
-      deallocate(wsumSTD, phisumSTD)
-      deallocate(tnondef)
-      deallocate(qsumSTD, rhsumSTD)
-      deallocate(uvsumSTD)
-      deallocate(vqsumSTD)
-      deallocate(vTsumSTD)
-      deallocate(wqsumSTD)
-      deallocate(vphisumSTD)
-      deallocate(wTsumSTD)
-      deallocate(u2sumSTD)
-      deallocate(v2sumSTD)
-      deallocate(T2sumSTD)
-      deallocate(O3sumSTD)
-      deallocate(O3daysumSTD)
+      DEALLOCATE(tsumSTD)
+      DEALLOCATE(usumSTD, vsumSTD)
+      DEALLOCATE(wsumSTD, phisumSTD)
+      DEALLOCATE(tnondef)
+      DEALLOCATE(qsumSTD, rhsumSTD)
+      DEALLOCATE(uvsumSTD)
+      DEALLOCATE(vqsumSTD)
+      DEALLOCATE(vTsumSTD)
+      DEALLOCATE(wqsumSTD)
+      DEALLOCATE(vphisumSTD)
+      DEALLOCATE(wTsumSTD)
+      DEALLOCATE(u2sumSTD)
+      DEALLOCATE(v2sumSTD)
+      DEALLOCATE(T2sumSTD)
+      DEALLOCATE(O3sumSTD)
+      DEALLOCATE(O3daysumSTD)
 !IM beg
-      deallocate(wlevSTD,ulevSTD,vlevSTD,tlevSTD,qlevSTD,rhlevSTD,philevSTD)
-      deallocate(uvSTD,vqSTD,vTSTD,wqSTD,vphiSTD,wTSTD,u2STD,v2STD,T2STD,O3STD,O3daySTD)
+      DEALLOCATE(wlevSTD,ulevSTD,vlevSTD,tlevSTD,qlevSTD,rhlevSTD,philevSTD)
+      DEALLOCATE(uvSTD,vqSTD,vTSTD,wqSTD,vphiSTD,wTSTD,u2STD,v2STD,T2STD,O3STD,O3daySTD)
 !IM end
-      deallocate(seed_old)
-      deallocate(zuthe, zvthe)
-      deallocate(alb_neig)
-      deallocate(ema_cbmf)
-      deallocate(ema_pcb, ema_pct)
-      deallocate(Mipsh, Ma, qcondc)
-      deallocate(wd, sigd)
-      deallocate(cin, ALE, ALP)
-      deallocate(ftd, fqd)
-      deallocate(Ale_bl, Alp_bl)
-      deallocate(ale_wake)
-      deallocate(ale_bl_stat)
-      deallocate(lalim_conv, wght_th)
-      deallocate(wake_deltat, wake_deltaq)
-      deallocate(wake_s, awake_dens, wake_dens)
-      deallocate(wake_Cstar, wake_pe, wake_fip)
+      DEALLOCATE(seed_old)
+      DEALLOCATE(zuthe, zvthe)
+      DEALLOCATE(alb_neig)
+      DEALLOCATE(ema_cbmf)
+      DEALLOCATE(ema_pcb, ema_pct)
+      DEALLOCATE(Mipsh, Ma, qcondc)
+      DEALLOCATE(wd, sigd)
+      DEALLOCATE(cin, ALE, ALP)
+      DEALLOCATE(ftd, fqd)
+      DEALLOCATE(Ale_bl, Alp_bl)
+      DEALLOCATE(ale_wake)
+      DEALLOCATE(ale_bl_stat)
+      DEALLOCATE(lalim_conv, wght_th)
+      DEALLOCATE(wake_deltat, wake_deltaq)
+      DEALLOCATE(wake_s, awake_dens, wake_dens)
+      DEALLOCATE(wake_Cstar, wake_pe, wake_fip)
 !jyg<
-      deallocate(wake_delta_pbl_TKE)
+      DEALLOCATE(wake_delta_pbl_TKE)
 !>jyg
-      deallocate(pfrac_impa, pfrac_nucl)
-      deallocate(pfrac_1nucl)
-      deallocate(total_rain, nday_rain)
-      deallocate(paire_ter)
-      deallocate(albsol1, albsol2)
+      DEALLOCATE(pfrac_impa, pfrac_nucl)
+      DEALLOCATE(pfrac_1nucl)
+      DEALLOCATE(total_rain, nday_rain)
+      DEALLOCATE(paire_ter)
+      DEALLOCATE(albsol1, albsol2)
 !albedo SB >>>
-      deallocate(albsol_dir,albsol_dif,falb_dir,falb_dif,chl_con)
+      DEALLOCATE(albsol_dir,albsol_dif,falb_dir,falb_dif,chl_con)
 !albedo SB <<<
-      deallocate(wo)
-      deallocate(clwcon0,rnebcon0)
-      deallocate(heat, heat0) 
-      deallocate(cool, cool0)
-      deallocate(heat_volc, cool_volc) 
-      deallocate(topsw, toplw)
-      deallocate(sollwdown, sollwdownclr)
-      deallocate(gustiness)
-      deallocate(toplwdown, toplwdownclr)
-      deallocate(topsw0,toplw0,solsw0,sollw0)
-      deallocate(albpla)
+      DEALLOCATE(wo)
+      DEALLOCATE(clwcon0,rnebcon0)
+      DEALLOCATE(heat, heat0) 
+      DEALLOCATE(cool, cool0)
+      DEALLOCATE(heat_volc, cool_volc) 
+      DEALLOCATE(topsw, toplw)
+      DEALLOCATE(sollwdown, sollwdownclr)
+      DEALLOCATE(gustiness)
+      DEALLOCATE(toplwdown, toplwdownclr)
+      DEALLOCATE(topsw0,toplw0,solsw0,sollw0)
+      DEALLOCATE(albpla)
 !IM ajout variables CFMIP2/CMIP5
-      deallocate(heatp, coolp)
-      deallocate(heat0p, cool0p)
-      deallocate(radsolp, topswp, toplwp)
-      deallocate(albplap)
-      deallocate(solswp, sollwp)
-      deallocate(sollwdownp)
-      deallocate(topsw0p,toplw0p)
-      deallocate(solsw0p,sollw0p)
-      deallocate(lwdnc0p, lwdn0p, lwdnp)
-      deallocate(lwupc0p, lwup0p, lwupp)
-      deallocate(swdnc0p, swdn0p, swdnp)
-      deallocate(swupc0p, swup0p, swupp)
-      deallocate(cape)
-      deallocate(pbase,bbase)
-      deallocate(zqasc)
-      deallocate(ibas_con, itop_con)
-      deallocate(rain_con, snow_con)
-      deallocate(rlonPOS)
-      deallocate(newsst)
-      deallocate(ustar,u10m, v10m,wstar)
-      deallocate(topswad, solswad)
-      deallocate(topswai, solswai)
-      deallocate(tau_aero,piz_aero,cg_aero)
-      deallocate(tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
-      deallocate(tau_aero_lw_rrtm,piz_aero_lw_rrtm,cg_aero_lw_rrtm)
-      deallocate(ccm)
-      if (ok_gwd_rando) deallocate(du_gwd_rando)
-      if (.not. ok_hines .and. ok_gwd_rando) deallocate(du_gwd_front)
+      DEALLOCATE(heatp, coolp)
+      DEALLOCATE(heat0p, cool0p)
+      DEALLOCATE(radsolp, topswp, toplwp)
+      DEALLOCATE(albplap)
+      DEALLOCATE(solswp, solswfdiffp, sollwp)
+      DEALLOCATE(sollwdownp)
+      DEALLOCATE(topsw0p,toplw0p)
+      DEALLOCATE(solsw0p,sollw0p)
+      DEALLOCATE(lwdnc0p, lwdn0p, lwdnp)
+      DEALLOCATE(lwupc0p, lwup0p, lwupp)
+      DEALLOCATE(swdnc0p, swdn0p, swdnp)
+      DEALLOCATE(swupc0p, swup0p, swupp)
+      DEALLOCATE(cape)
+      DEALLOCATE(pbase,bbase)
+      DEALLOCATE(zqasc)
+      DEALLOCATE(ibas_con, itop_con)
+      DEALLOCATE(rain_con, snow_con)
+      DEALLOCATE(rlonPOS)
+      DEALLOCATE(newsst)
+      DEALLOCATE(ustar,u10m, v10m,wstar)
+      DEALLOCATE(topswad, solswad)
+      DEALLOCATE(topswai, solswai)
+      DEALLOCATE(tau_aero,piz_aero,cg_aero)
+      DEALLOCATE(tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
+      DEALLOCATE(tau_aero_lw_rrtm,piz_aero_lw_rrtm,cg_aero_lw_rrtm)
+      DEALLOCATE(ccm)
+      if (ok_gwd_rando) DEALLOCATE(du_gwd_rando)
+      if (.not. ok_hines .and. ok_gwd_rando) DEALLOCATE(du_gwd_front)
        
 !!! nrlmd le 10/04/2012
-      deallocate(ale_bl_trig)
+      DEALLOCATE(ale_bl_trig)
 !!! fin nrlmd le 10/04/2012
       is_initialized=.FALSE.
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3756)
@@ -2514,5 +2514,5 @@
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
             zsig,      sollwdown, pphi,    cldt,      &
-            rain_fall, snow_fall, solsw,   sollw,     &
+            rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
             gustiness,                                &
             t_seri,    q_seri,    u_seri,  v_seri,    &
@@ -4143,5 +4143,5 @@
                heat,heat0,cool,cool0,albpla, &
                heat_volc,cool_volc, &
-               topsw,toplw,solsw,sollw, &
+               topsw,toplw,solsw,solswfdiff,sollw, &
                sollwdown, &
                topsw0,toplw0,solsw0,sollw0, &
@@ -4230,5 +4230,5 @@
                      heatp,heat0p,coolp,cool0p,albplap, &
                      heat_volc,cool_volc, &
-                     topswp,toplwp,solswp,sollwp, &
+                     topswp,toplwp,solswp,solswfdiffp,sollwp, &
                      sollwdownp, &
                      topsw0p,toplw0p,solsw0p,sollw0p, &
Index: LMDZ6/trunk/libf/phylmd/radlwsw_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/radlwsw_m.F90	(revision 3748)
+++ LMDZ6/trunk/libf/phylmd/radlwsw_m.F90	(revision 3756)
@@ -26,5 +26,5 @@
    heat,heat0,cool,cool0,albpla,&
    heat_volc, cool_volc,&
-   topsw,toplw,solsw,sollw,&
+   topsw,toplw,solsw,solswfdiff,sollw,&
    sollwdown,&
    topsw0,toplw0,solsw0,sollw0,&
@@ -117,4 +117,5 @@
   ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
   ! solsw----output-R- flux solaire net a la surface
+  ! solswfdiff----output-R- fraction de rayonnement diffus pour le flux solaire descendant a la surface
   ! sollw----output-R- ray. IR montant a la surface
   ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
@@ -188,5 +189,5 @@
   REAL,    INTENT(in)  :: tsol(KLON)
   REAL,    INTENT(in) :: alb_dir(KLON,NSW),alb_dif(KLON,NSW)
-  real, intent(in) :: SFRWL(6)
+  REAL,    INTENT(in) :: SFRWL(6)
 !albedo SB <<<
   REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV)
@@ -235,5 +236,5 @@
   REAL,    INTENT(out) :: heat_volc(KLON,KLEV), cool_volc(KLON,KLEV) !NL
   REAL,    INTENT(out) :: topsw(KLON), toplw(KLON)
-  REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON)
+  REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON), solswfdiff(KLON)
   REAL,    INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
   REAL,    INTENT(out) :: sollwdown(KLON)
@@ -286,5 +287,5 @@
   REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
 
-  real(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
+  REAL(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
   ! "POZON(:, :, 1)" is for the average day-night field, 
   ! "POZON(:, :, 2)" is for daylight time.
@@ -302,5 +303,5 @@
   REAL(KIND=8) zheat_volc(kdlon,kflev), zcool_volc(kdlon,kflev) !NL
   REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
-  REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
+  REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon), zsolswfdiff(kdlon)
   REAL(KIND=8) zsollwdown(kdlon)
   REAL(KIND=8) ztopsw0(kdlon), ztoplw0(kdlon)
@@ -329,74 +330,76 @@
 !MPL input supplementaires pour RECMWFL
 ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
-      REAL(KIND=8) GEMU(klon)
+  REAL(KIND=8) GEMU(klon)
 !MPL input RECMWFL: 
 ! Tableaux aux niveaux inverses pour respecter convention Arpege
-      REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
-      REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
+  REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
+  REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
 !--OB
-      REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
-      REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
+  REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
+  REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
 !--end OB
-      REAL(KIND=8) paprs_i(klon,klev+1)
-      REAL(KIND=8) pplay_i(klon,klev)
-      REAL(KIND=8) cldfra_i(klon,klev)
-      REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone 
+  REAL(KIND=8) paprs_i(klon,klev+1)
+  REAL(KIND=8) pplay_i(klon,klev)
+  REAL(KIND=8) cldfra_i(klon,klev)
+  REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone 
   ! "POZON(:, :, 1)" is for the average day-night field, 
   ! "POZON(:, :, 2)" is for daylight time.
 !!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6      
-      REAL(KIND=8) PAER_i(kdlon,kflev,6)
-      REAL(KIND=8) PDP_i(klon,klev)
-      REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
-      REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
+  REAL(KIND=8) PAER_i(kdlon,kflev,6)
+  REAL(KIND=8) PDP_i(klon,klev)
+  REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
+  REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
 !MPL output RECMWFL:
-      REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)       
-      REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)      
-      REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)    
-      REAL(KIND=8) ZTH   (klon,klev+1),ZTH_i   (klon,klev+1)   
-      REAL(KIND=8) ZCTRSO(klon,2)       
-      REAL(KIND=8) ZCEMTR(klon,2)     
-      REAL(KIND=8) ZTRSOD(klon)        
-      REAL(KIND=8) ZLWFC (klon,2)      
-      REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)    
-      REAL(KIND=8) ZSWFC (klon,2)      
-      REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
-      REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
-      REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
-      REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
-      REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
-      REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
-      REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
-      REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
+  REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)       
+  REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)      
+  REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)    
+  REAL(KIND=8) ZTH   (klon,klev+1),ZTH_i   (klon,klev+1)   
+  REAL(KIND=8) ZCTRSO(klon,2)       
+  REAL(KIND=8) ZCEMTR(klon,2)     
+  REAL(KIND=8) ZTRSOD(klon)        
+  REAL(KIND=8) ZLWFC (klon,2)      
+  REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)    
+  REAL(KIND=8) ZSWFC (klon,2)      
+  REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
+  REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
+  REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
+  REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
+  REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
+  REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
+  REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
+  REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
 #ifdef CPP_RRTM
-      REAL(KIND=8) PTAU_LW_TOT(klon,klev,NLW)
-      REAL(KIND=8) PTAU_LW_NAT(klon,klev,NLW)
+  REAL(KIND=8) PTAU_LW_TOT(klon,klev,NLW)
+  REAL(KIND=8) PTAU_LW_NAT(klon,klev,NLW)
 #endif
-      REAL(KIND=8) PSFSWDIR(klon,NSW)
-      REAL(KIND=8) PSFSWDIF(klon,NSW)
-      REAL(KIND=8) PFSDNN(klon)
-      REAL(KIND=8) PFSDNV(klon)
+  REAL(KIND=8) PSFSWDIR(klon,NSW)
+  REAL(KIND=8) PSFSWDIF(klon,NSW)
+  REAL(KIND=8) PFSDNN(klon)
+  REAL(KIND=8) PFSDNV(klon)
 !MPL On ne redefinit pas les tableaux ZFLUX,ZFLUC,
 !MPL ZFSDWN,ZFCDWN,ZFSUP,ZFCUP car ils existent deja
 !MPL sous les noms de ZFLDN,ZFLDN0,ZFLUP,ZFLUP0,
 !MPL ZFSDN,ZFSDN0,ZFSUP,ZFSUP0
-      REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
-      REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
-      REAL(KIND=8) ZFSDWN_i (klon,klev+1)
-      REAL(KIND=8) ZFCDWN_i (klon,klev+1)
-      REAL(KIND=8) ZFCCDWN_i (klon,klev+1)
-      REAL(KIND=8) ZFSUP_i (klon,klev+1)
-      REAL(KIND=8) ZFCUP_i (klon,klev+1)
-      REAL(KIND=8) ZFCCUP_i (klon,klev+1)
-      REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
-      REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
+  REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
+  REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
+  REAL(KIND=8) ZFSDWN_i (klon,klev+1)
+  REAL(KIND=8) ZFCDWN_i (klon,klev+1)
+  REAL(KIND=8) ZFCCDWN_i (klon,klev+1)
+  REAL(KIND=8) ZFSUP_i (klon,klev+1)
+  REAL(KIND=8) ZFCUP_i (klon,klev+1)
+  REAL(KIND=8) ZFCCUP_i (klon,klev+1)
+  REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
+  REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
 ! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
 !      REAL(KIND=8) RSUN(3,2)
 !      REAL(KIND=8) SUN(3)
 !      REAL(KIND=8) SUN_FRACT(2)
-  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
+  REAL, PARAMETER:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
   CHARACTER (LEN=80) :: abort_message
   CHARACTER (LEN=80) :: modname='radlwsw_m'
 
-  call assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
+  REAL zdir, zdif
+
+  CALL assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
   ! initialisation
   ist=1
@@ -414,13 +417,12 @@
   zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
 
-
-   ZTOPSWADAERO(:)  = 0. !ym missing init
-   ZSOLSWADAERO(:)  = 0. !ym missing init
-   ZTOPSWAD0AERO(:) = 0. !ym missing init
-   ZSOLSWAD0AERO(:) = 0. !ym missing init
-   ZTOPSWAIAERO(:)  = 0. !ym missing init 
-   ZSOLSWAIAERO(:)  = 0. !ym missing init  
-   ZTOPSWCF_AERO(:,:)= 0.!ym missing init  
-   ZSOLSWCF_AERO(:,:) =0. !ym missing init  
+  ZTOPSWADAERO(:)  = 0. !ym missing init
+  ZSOLSWADAERO(:)  = 0. !ym missing init
+  ZTOPSWAD0AERO(:) = 0. !ym missing init
+  ZSOLSWAD0AERO(:) = 0. !ym missing init
+  ZTOPSWAIAERO(:)  = 0. !ym missing init 
+  ZSOLSWAIAERO(:)  = 0. !ym missing init  
+  ZTOPSWCF_AERO(:,:)= 0.!ym missing init  
+  ZSOLSWCF_AERO(:,:) =0. !ym missing init  
 
   !
@@ -454,9 +456,9 @@
 #ifdef REPROBUS
     IF (iflag_rrtm==0) THEN
-      if(ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
+      IF (ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
       print*,'Constante solaire: ',PSCT*zdist*zdist
-    END IF
+    ENDIF
 #endif
-  END IF
+  ENDIF
 
   DO j = 1, nb_gr
@@ -541,6 +543,5 @@
        CALL RAD_INTERACTIF(POZON,iof)
 #endif
-    END IF
-
+    ENDIF
     !
     DO k = 1, kflev+1
@@ -568,9 +569,9 @@
       ENDDO
     ENDDO
-
 !
 !===== iflag_rrtm ================================================
 !      
     IF (iflag_rrtm == 0) THEN       !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
+!
 !--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------              
       DO k = 1, kflev+1
@@ -651,4 +652,9 @@
       zsolswaiaero(i)=0.
       ENDDO
+
+      !--fraction of diffuse radiation in surface SW downward radiation
+      !--not computed with old radiation scheme
+      zsolswfdiff(:) = -999.999
+
 !     print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
        ! daylight ozone, if we have it, for short wave
@@ -676,5 +682,4 @@
        DO i=1,kdlon
        DO k=1,kflev+1
-!        print *,'iof i k klon klev=',iof,i,k,klon,klev
          lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
          lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
@@ -687,10 +692,5 @@
        ENDDO  
        ENDDO  
-!          print*,'SW_AR4 ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
-!          print*,'SW_AR4 swdn0  1 , klev:',swdn0(1:klon,1),swdn0(1:klon,klev)
-!          print*,'SW_AR4 ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
-!          print*,'SW_AR4 swup0  1 , klev:',swup0(1:klon,1),swup0(1:klon,klev)
-!          print*,'SW_AR4 ZFSDN  1 , klev:',ZFSDN(1:klon,1) ,ZFSDN(1:klon,klev)
-!          print*,'SW_AR4 ZFSUP  1 , klev:',ZFSUP(1:klon,1) ,ZFSUP(1:klon,klev)
+!
     ELSE  
 #ifdef CPP_RRTM
@@ -700,22 +700,22 @@
       DO k = 1, kflev+1
       DO i = 1, kdlon
-      ZEMTD_i(i,k)=0.
-      ZEMTU_i(i,k)=0.
-      ZTRSO_i(i,k)=0.
-      ZTH_i(i,k)=0.
-      ZLWFT_i(i,k)=0.
-      ZSWFT_i(i,k)=0.
-      ZFLUX_i(i,1,k)=0.
-      ZFLUX_i(i,2,k)=0.
-      ZFLUC_i(i,1,k)=0.
-      ZFLUC_i(i,2,k)=0.
-      ZFSDWN_i(i,k)=0.
-      ZFCDWN_i(i,k)=0.
-      ZFCCDWN_i(i,k)=0.
-      ZFSUP_i(i,k)=0.
-      ZFCUP_i(i,k)=0.
-      ZFCCUP_i(i,k)=0.
-      ZFLCCDWN_i(i,k)=0.
-      ZFLCCUP_i(i,k)=0.
+        ZEMTD_i(i,k)=0.
+        ZEMTU_i(i,k)=0.
+        ZTRSO_i(i,k)=0.
+        ZTH_i(i,k)=0.
+        ZLWFT_i(i,k)=0.
+        ZSWFT_i(i,k)=0.
+        ZFLUX_i(i,1,k)=0.
+        ZFLUX_i(i,2,k)=0.
+        ZFLUC_i(i,1,k)=0.
+        ZFLUC_i(i,2,k)=0.
+        ZFSDWN_i(i,k)=0.
+        ZFCDWN_i(i,k)=0.
+        ZFCCDWN_i(i,k)=0.
+        ZFSUP_i(i,k)=0.
+        ZFCUP_i(i,k)=0.
+        ZFCCUP_i(i,k)=0.
+        ZFLCCDWN_i(i,k)=0.
+        ZFLCCUP_i(i,k)=0.
       ENDDO
       ENDDO
@@ -771,6 +771,6 @@
       PFSDNV(i)=0.
       DO kk = 1, NSW
-      PSFSWDIR(i,kk)=0.
-      PSFSWDIF(i,kk)=0.
+        PSFSWDIR(i,kk)=0.
+        PSFSWDIF(i,kk)=0.
       ENDDO
       ENDDO
@@ -779,5 +779,5 @@
 ! On met les donnees dans l'ordre des niveaux arpege
          paprs_i(:,1)=paprs(:,klev+1)
-         do k=1,klev
+         DO k=1,klev
             paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
             pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
@@ -794,14 +794,14 @@
             ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
             ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
-         enddo
-         do k=1,kflev
+         ENDDO
+         DO k=1,kflev
            POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
 !!!            POZON_i(1:klon,k)=POZON(1:klon,k)	    !!! on laisse 1=sol et klev=top 
 !          print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
 !!!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6      
-            do i=1,6
+            DO i=1,6
             PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
-            enddo
-         enddo
+            ENDDO
+         ENDDO
 !       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
 
@@ -832,5 +832,5 @@
 !    s   'RECMWF ')
 !
-      if(lldebug) then
+      IF (lldebug) THEN
         CALL writefield_phy('paprs_i',paprs_i,klev+1)
         CALL writefield_phy('pplay_i',pplay_i,klev)
@@ -846,5 +846,5 @@
         CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
         CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
-      endif
+      ENDIF
 
 ! Nouvel appel a RECMWF (celui du cy32t0)
@@ -876,5 +876,5 @@
             
 !        print *,'RADLWSW: apres RECMWF'
-      if(lldebug) then
+      IF (lldebug) THEN
         CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
         CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
@@ -901,5 +901,5 @@
         CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
         CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
-      endif
+      ENDIF
 ! --------- output RECMWFL
 !  ZEMTD        (KPROMA,KLEV+1)  ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
@@ -952,5 +952,5 @@
          ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
          ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
-         IF(ok_volcan) THEN
+         IF (ok_volcan) THEN
             ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
          ENDIF
@@ -992,4 +992,16 @@
 ! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
 ! LW_LMDAR4 et SW_LMDAR4
+
+      !--fraction of diffuse radiation in surface SW downward radiation
+      DO i = 1, kdlon
+       IF (fract(i).GT.0.0) THEN
+         zdir=SUM(PSFSWDIR(i,:))
+         zdif=SUM(PSFSWDIF(i,:))
+         zsolswfdiff(i) = zdif/(zdir+zdif)
+       ELSE  !--night
+         zsolswfdiff(i) = 1.0
+       ENDIF
+      ENDDO 
+!
       DO i = 1, kdlon
          zsolsw(i)    = ZSWFT(i,1)
@@ -1009,5 +1021,5 @@
          ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
 !         
-           IF (fract(i) == 0.) THEN
+         IF (fract(i) == 0.) THEN
 !!!!! A REVOIR MPL (20090630) ca n a pas de sens quand fract=0
 ! pas plus que dans le sw_AR4
@@ -1030,11 +1042,11 @@
 ! ZLWFT(klon,k),ZSWFT
 
-      do k=1,kflev
-         do i=1,kdlon
+      DO k=1,kflev
+         DO i=1,kdlon
            zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
            zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
            zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
            zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
-           IF(ok_volcan) THEN
+           IF (ok_volcan) THEN
               zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
               zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
@@ -1043,6 +1055,6 @@
 !	   ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
 !	   ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)	   
-         enddo
-      enddo
+         ENDDO
+      ENDDO
 #else
     abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
@@ -1056,4 +1068,5 @@
       toplw(iof+i) = ztoplw(i)
       solsw(iof+i) = zsolsw(i)
+      solswfdiff(iof+i) = zsolswfdiff(i)
       sollw(iof+i) = zsollw(i)
       sollwdown(iof+i) = zsollwdown(i)
