Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clcdrag.F90
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clcdrag.F90	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clcdrag.F90	(revision 467)
@@ -40,5 +40,6 @@
 !
 ! Quelques constantes et options:
-      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
+!!$PB      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
+      REAL, PARAMETER :: ckap=0.40, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
 !
 ! Variables locales :
@@ -68,5 +69,6 @@
         zcdn(i) = (ckap/log(1.+zgeop(i)/(RG*rugos(i))))**2
 !
-        IF (zri(i) .ge. 0.) THEN      ! situation stable
+!!$        IF (zri(i) .ge. 0.) THEN      ! situation stable
+        IF (zri(i) .gt. 0.) THEN      ! situation stable
           zri(i) = min(20.,zri(i))
           IF (.NOT.zxli) THEN
@@ -75,5 +77,6 @@
             zcfm1(i) = zcdn(i) * FRIV
             FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
-            zcfh1(i) = zcdn(i) * FRIH
+!!$  PB          zcfh1(i) = zcdn(i) * FRIH
+            zcfh1(i) = 0.8 * zcdn(i) * FRIH
             pcfm(i) = zcfm1(i)
             pcfh(i) = zcfh1(i)
@@ -87,5 +90,6 @@
                  *(1.0+zgeop(i)/(RG*rugos(i)))))
             zcfm2(i) = zcdn(i)*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
-            zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
+!!$PB            zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
+            zcfh2(i) = 0.8 * zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
             pcfm(i) = zcfm2(i)
             pcfh(i) = zcfh2(i)
@@ -94,6 +98,6 @@
             pcfh(i) = zcdn(i)* fins(zri(i))
           ENDIF
-          zcr = (0.0016/(zcdn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
-          IF(nsrf.EQ.is_oce) pcfh(i) = zcdn(i)*(1.0+zcr**1.25)**(1./1.25)
+  	          zcr = (0.0016/(zcdn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
+          IF(nsrf.EQ.is_oce) pcfh(i) =0.8* zcdn(i)*(1.0+zcr**1.25)**(1./1.25)
         ENDIF
       END DO
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.h	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.h	(revision 467)
@@ -6,6 +6,13 @@
        REAL co2_ppm, solaire
        REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12  
+       REAL*8 CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt
+cIM simulateur ISCCP 
+       INTEGER top_height, overlap
+cIM seuils cdrm, cdrh
+       REAL cdmmax, cdhmax
 
        COMMON/clesphys/cycle_diurne, soil_model, new_oliq,
      ,     ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad, iflag_con
      ,     , co2_ppm, solaire, RCO2, RCH4, RN2O, RCFC11, RCFC12
+     ,     , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt
+     ,     , top_height, overlap, cdmmax, cdhmax
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.inc
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.inc	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.inc	(revision 467)
@@ -6,6 +6,11 @@
        REAL :: co2_ppm, solaire
        DOUBLE PRECISION :: RCO2, RCH4, RN2O, RCFC11, RCFC12  
+       DOUBLE PRECISION :: CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt
+       INTEGER :: top_height, overlap
+       REAL :: cdmmax, cdhmax
 
        COMMON/clesphys/cycle_diurne, soil_model, new_oliq, &
      &     ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad, iflag_con &
-     &     , co2_ppm, solaire, RCO2, RCH4, RN2O, RCFC11, RCFC12
+     &     , co2_ppm, solaire, RCO2, RCH4, RN2O, RCFC11, RCFC12 &
+     &     , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt &
+     &     , top_height, overlap, cdmmax, cdhmax
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 467)
@@ -7,5 +7,5 @@
      .                  jour, rmu0,
      .                  ok_veget, ocean, npas, nexca, ts,
-     .                  soil_model,ftsoil,qsol,
+     .                  soil_model,cdmmax, cdhmax, ftsoil,qsol,
      .                  paprs,pplay,radsol,snow,qsurf,evap,albe,alblw,
      .                  fluxlat,
@@ -140,4 +140,6 @@
 c$$$ PB ajout pour soil
       LOGICAL soil_model
+cIM ajout seuils cdrm, cdrh
+      REAL cdmmax, cdhmax
       REAL ftsoil(klon,nsoilmx,nbsrf)
       REAL ytsoil(klon,nsoilmx)
@@ -481,5 +483,14 @@
       ENDDO
 c
-c
+cIM cf JLD : on seuille ycoefm et ycoefh
+      if (nsrf.eq.is_oce) then
+         do j=1,knon
+c           ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)
+            ycoefm(j,1)=min(ycoefm(j,1),cdmmax)
+c           ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)
+            ycoefh(j,1)=min(ycoefh(j,1),cdhmax)
+         enddo
+      endif
+
 c calculer la diffusion des vitesses "u" et "v"
       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
@@ -493,7 +504,8 @@
 
 c FH modif sur le cdrag temperature
-      do i=1,knon
-         ycoefh(i,1)=ycoefm(i,1)*0.8
-      enddo
+c$$$PB : déplace dans clcdrag
+c$$$      do i=1,knon
+c$$$         ycoefh(i,1)=ycoefm(i,1)*0.8
+c$$$      enddo
 
 c calculer la diffusion de "q" et de "h"
@@ -520,5 +532,6 @@
       IF (nsrf.EQ.is_oce) THEN
       DO j = 1, knon
-         yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
+         yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG 
+     $      +  0.11*14e-6 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))
          yrugm(j) = MAX(1.5e-05,yrugm(j))
       ENDDO
@@ -1286,4 +1299,7 @@
       REAL t_coup
       PARAMETER (t_coup=273.15)
+cIM
+      LOGICAL check
+      PARAMETER (check=.false.)
 c
 c contre-gradient pour la chaleur sensible: Kelvin/metre
@@ -1388,5 +1404,7 @@
       ENDDO
 
-c$$$      PRINT*,' isommet=',isommet,' knon=',knon
+      IF (check) THEN
+      PRINT*,' isommet=',isommet,' knon=',knon
+      ENDIF
 
       DO k = 2, isommet
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 467)
@@ -168,10 +168,13 @@
 !
 !valeur AMIP II
-  RCH4 = 1.65E-06* 16.043/28.97
+!OK  RCH4 = 1.65E-06* 16.043/28.97
 ! RCH4 = 9.137366240938903E-07
 !
 !ancienne valeur
 ! RCH4 = 1.72E-06* 16.043/28.97
-  call getin('RCH4', RCH4)
+!OK call getin('RCH4', RCH4)
+  CH4_ppb = 1650.
+  call getin('CH4_ppb', CH4_ppb)
+  RCH4 = CH4_ppb * 1.0E-09 * 16.043/28.97
 !!
 !Config Key  = RN2O
@@ -182,10 +185,13 @@
 !
 !valeur AMIP II
-  RN2O = 306.E-09* 44.013/28.97
+!OK  RN2O = 306.E-09* 44.013/28.97
 ! RN2O = 4.648939592682085E-07
 !
 !ancienne valeur
 ! RN2O = 310.E-09* 44.013/28.97
-  call getin('RN2O', RN2O)
+!OK  call getin('RN2O', RN2O)
+  N2O_ppb=306.
+  call getin('N2O_ppb', N2O_ppb)
+  RN2O = N2O_ppb * 1.0E-09 * 44.013/28.97
 !!
 !Config Key  = RCFC11
@@ -195,7 +201,10 @@
 !               
 !
-  RCFC11 = 280.E-12* 137.3686/28.97
+!OK RCFC11 = 280.E-12* 137.3686/28.97
+  CFC11_ppt = 280.
+  call getin('CFC11_ppt',CFC11_ppt)
+  RCFC11=CFC11_ppt* 1.0E-12 * 137.3686/28.97
 ! RCFC11 = 1.327690990680013E-09
-  call getin('RCFC11', RCFC11)
+!OK call getin('RCFC11', RCFC11)
 !!
 !Config Key  = RCFC12
@@ -205,7 +214,10 @@
 !               
 !
-  RCFC12 = 484.E-12* 120.9140/28.97
+!OK RCFC12 = 484.E-12* 120.9140/28.97
+  CFC12_ppt = 484.
+  call getin('CFC12_ppt',CFC12_ppt)
+  RCFC12 = CFC12_ppt * 1.0E-12 * 120.9140/28.97
 ! RCFC12 = 2.020102726958923E-09
-  call getin('RCFC12', RCFC12)
+!OK call getin('RCFC12', RCFC12)
 !!
 !! Constante solaire & Parametres orbitaux & taux gaz effet de serre END
@@ -386,4 +398,42 @@
   rad_chau2 = 9.0
   call getin('rad_chau2',rad_chau2)
+
+!
+!Config Key  = top_height
+!Config Desc =
+!Config Def  = 3
+!Config Help =
+!
+  top_height = 3
+  call getin('top_height',top_height)
+
+!
+!Config Key  = overlap
+!Config Desc =
+!Config Def  = 3
+!Config Help =
+!
+  overlap = 3
+  call getin('overlap',overlap)
+
+
+!IM
+!
+!Config Key  = cdmmax
+!Config Desc =
+!Config Def  = 1.3E-3
+!Config Help =
+!
+  cdmmax = 1.3E-3
+  call getin('cdmmax',cdmmax)
+
+!
+!Config Key  = cdhmax
+!Config Desc =
+!Config Def  = 1.1E-3
+!Config Help =
+!
+  cdhmax = 1.1E-3
+  call getin('cdhmax',cdhmax)
 
 !
@@ -415,8 +465,8 @@
   write(numout,*)' co2_ppm =',co2_ppm
   write(numout,*)' RCO2 = ',RCO2
-  write(numout,*)' RCH4 = ',RCH4
-  write(numout,*)' RN2O =  ',RN2O
-  write(numout,*)' RCFC11 =  ',RCFC11
-  write(numout,*)' RCFC12 =  ',RCFC12
+  write(numout,*)' CH4_ppb =',CH4_ppb,' RCH4 = ',RCH4
+  write(numout,*)' N2O_ppb =',N2O_ppb,' RN2O =  ',RN2O
+  write(numout,*)' CFC11_ppt=',CFC11_ppt,' RCFC11 =  ',RCFC11
+  write(numout,*)' CFC12_ppt=',CFC12_ppt,' RCFC12 =  ',RCFC12
 !IM constantes physiques END
   write(numout,*)' epmax = ', epmax
@@ -438,4 +488,6 @@
   write(numout,*)' ratqsbas = ',ratqsbas 
   write(numout,*)' ratqshaut = ',ratqshaut 
+  write(numout,*)' top_height = ',top_height 
+  write(numout,*)' overlap = ',overlap 
 
   return
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histday.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histday.h	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histday.h	(revision 467)
@@ -170,4 +170,28 @@
      .                32, "ave(X)", zsto,zout)
 c
+         CALL histdef(nid_day, "SWupTOAclr",
+     .                "SWup clear sky at TOA","W/m2",
+     .                iim,jjmp1,nhori, 1,1,1,-99,
+     .                32, "ave(X)", zsto,zout)
+
+         CALL histdef(nid_day, "SWupSFCclr",
+     .                "SWup clear sky at surface","W/m2",
+     .                iim,jjmp1,nhori, 1,1,1,-99,
+     .                32, "ave(X)", zsto,zout)
+
+         CALL histdef(nid_day, "SWdnTOAclr",
+     .                "SWdn clear sky at TOA","W/m2",
+     .                iim,jjmp1,nhori, 1,1,1,-99,
+     .                32, "ave(X)", zsto,zout)
+
+         CALL histdef(nid_day, "SWdnSFCclr",
+     .                "SWdn clear sky at surface","W/m2",
+     .                iim,jjmp1,nhori, 1,1,1,-99,
+     .                32, "ave(X)", zsto,zout)
+cccIM   
+         CALL histdef(nid_day, "prw", "Precipitable water", "kg/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "ave(X)", zsto,zout)
+c
 c  Champs dynamiques sur niveaux de pression
 
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histhf.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histhf.h	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histhf.h	(revision 467)
@@ -83,8 +83,25 @@
      .                "inst(X)", zsto,zout)
 
-         CALL histdef(nid_hf, "phi500", "Geopotentiel à 500mb", "m2/s2",
+         CALL histdef(nid_hf, "phi500", "Geopotentiel a 500mb", "m2/s2",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
      .                "inst(X)", zsto,zout)
 
+cIM cf FH
+         CALL histdef(nid_hf,"u1","Zonal wind at 1st layer", "m/s",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)", zsto,zout)
+
+         CALL histdef(nid_hf,"v1","Meridional wind at 1st layer",
+     .                "m/s",iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)", zsto,zout)
+ 
+         CALL histdef(nid_hf, "cdrm", " Momentum drag coef.", "-",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)", zsto,zout)
+
+         CALL histdef(nid_hf, "cdrh", "Heat drag coef.", "-",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)", zsto,zout)
+ 
 c
          CALL histend(nid_hf)
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h	(revision 467)
@@ -539,4 +539,9 @@
      .                iim,jjmp1,nhori, 1,1,1,-99,
      .                32, "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "prw", "Precipitable water", "kg/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
+     .                "ave(X)", zsto,zout)
+
 c Champs interpolles sur des niveaux de pression
 
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 467)
@@ -264,4 +264,8 @@
 !!$  tsurf_new = 0.
 
+!IM cf JLD
+  ffonte(1:knon)=0.
+  fqcalving(1:knon)=0.
+
   cal = 999999. ; beta = 999999. ; dif_grnd = 999999. ; capsol = 999999.
   alb_new = 999999. ; z0_new = 999999. ; alb_neig = 999999.
@@ -1096,7 +1100,9 @@
 !IM cf. JP ---
 
-   where(cdrag > 0.01) 
-     cdrag = 0.01
-   endwhere
+
+! PF et PASB
+!   where(cdrag > 0.01) 
+!     cdrag = 0.01
+!   endwhere
 !  write(*,*)'Cdrag = ',minval(cdrag),maxval(cdrag)
 
@@ -1116,5 +1122,5 @@
 
 !IM cf. JP +++
-    albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2.
+    albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
 !IM cf. JP ---
 
@@ -1122,5 +1128,5 @@
 
 !IM cf. JP +++
-  swdown_vrai(:) = swnet(:)/(1. - albedo_keep(:))
+  swdown_vrai(1:knon) = swnet(1:knon)/(1. - albedo_keep(1:knon))
 !IM cf. JP ---
 
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 467)
@@ -40,6 +40,5 @@
 #include "mpiclim.h"
 c
-#include "oasis.h"      
-                        ! contains the name of communication technique. Here
+#include "oasis.h"      ! contains the name of communication technique. Here
                         ! cchan=CLIM only is possible.
 c			! ctype=MPI2
@@ -555,19 +554,4 @@
       END
 
-      SUBROUTINE halte
-      print *, 'Attention dans oasis.F, halte est non defini'
-      RETURN
-      END
-
-      SUBROUTINE locread
-      print *, 'Attention dans oasis.F, locread est non defini'
-      RETURN
-      END
-
-      SUBROUTINE locwrite
-      print *, 'Attention dans oasis.F, locwrite est non defini'
-      RETURN
-      END
-
       SUBROUTINE pipe_model_define
       print*,'Attention dans oasis.F, pipe_model_define est non defini'
@@ -590,54 +574,4 @@
       END
 
-      SUBROUTINE clim_stepi
-      print *, 'Attention dans oasis.F, clim_stepi est non defini'
-      RETURN
-      END
-
-      SUBROUTINE clim_start
-      print *, 'Attention dans oasis.F, clim_start est non defini'
-      RETURN
-      END
-
-      SUBROUTINE clim_import
-      print *, 'Attention dans oasis.F, clim_import est non defini'
-      RETURN
-      END
-
-      SUBROUTINE clim_export
-      print *, 'Attention dans oasis.F, clim_export est non defini'
-      RETURN
-      END
-
-      SUBROUTINE clim_init
-      print *, 'Attention dans oasis.F, clim_init est non defini'
-      RETURN
-      END
-
-      SUBROUTINE clim_define
-      print *, 'Attention dans oasis.F, clim_define est non defini'
-      RETURN
-      END
-
-      SUBROUTINE clim_quit
-      print *, 'Attention dans oasis.F, clim_quit est non defini'
-      RETURN
-      END
-
-      SUBROUTINE svipc_write
-      print *, 'Attention dans oasis.F, svipc_write est non defini'
-      RETURN
-      END
-
-      SUBROUTINE svipc_close
-      print *, 'Attention dans oasis.F, svipc_close est non defini'
-      RETURN
-      END
-
-      SUBROUTINE svipc_read
-      print *, 'Attention dans oasis.F, svipc_read est non defini'
-      RETURN
-      END
-
       SUBROUTINE quitcpl
       print *, 'Attention dans oasis.F, quitcpl est non defini'
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F	(revision 467)
@@ -4,5 +4,5 @@
       SUBROUTINE phyetat0 (fichnom,dtime,co2_ppm_etat0,solaire_etat0,
      .            rlat,rlon, pctsrf, tsol,tsoil,deltat,qsurf,qsol,snow,
-     .           albe, evap, rain_fall, snow_fall, solsw, sollw,
+     .           albe, alblw, evap, rain_fall, snow_fall, solsw, sollw,
      .           fder,radsol,frugs,agesno,clesphy0,
      .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,tabcntr0,
@@ -34,4 +34,7 @@
       REAL snow(klon,nbsrf)
       REAL albe(klon,nbsrf)
+cIM BEG alblw
+      REAL alblw(klon,nbsrf)
+cIM END alblw
       REAL evap(klon,nbsrf)
       REAL radsol(klon)
@@ -645,4 +648,73 @@
 
 c
+cIM BEG alblw
+c Lecture de albedo au sol LW:
+c
+      ierr = NF_INQ_VARID (nid, "ALBLW", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <ALBLW> est absent'
+c        PRINT*, '          Mais je vais essayer de lire ALBLW**'
+         PRINT*, '          Mais je vais prendre ALBE**'
+         DO nsrf = 1, nbsrf
+           DO i = 1, klon
+             alblw(i,nsrf) = albe(i,nsrf)
+           ENDDO
+         ENDDO
+c          IF (nsrf.GT.99) THEN
+c            PRINT*, "Trop de sous-mailles"
+c            CALL abort
+c          ENDIF
+c          WRITE(str2,'(i2.2)') nsrf
+c          ierr = NF_INQ_VARID (nid, "ALBLW"//str2, nvarid)
+c           IF (ierr.NE.NF_NOERR) THEN
+c             PRINT*, "phyetat0: Le champ <ALBLW"//str2//"> est absent"
+c             CALL abort
+c          ENDIF
+c#ifdef NC_DOUBLE
+c           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alblw(1,nsrf))
+c#else
+c           ierr = NF_GET_VAR_REAL(nid, nvarid, alblw(1,nsrf))
+c#endif
+c          IF (ierr.NE.NF_NOERR) THEN
+c            PRINT*, "phyetat0: Lecture echouee pour <ALBLW"//str2//">"
+c            CALL abort
+c          ENDIF
+c          xmin = 1.0E+20
+c          xmax = -1.0E+20
+c          DO i = 1, klon
+c             xmin = MIN(alblw(i,nsrf),xmin)
+c             xmax = MAX(alblw(i,nsrf),xmax)
+c          ENDDO
+c          PRINT*,'Albedo du sol ALBLW**:', nsrf, xmin, xmax
+c        ENDDO
+      ELSE
+         PRINT*, 'phyetat0: Le champ <ALBLW> est present'
+         PRINT*, '          J ignore donc les autres ALBLW**'
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alblw(1,1))
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, alblw(1,1))
+#endif
+         IF (ierr.NE.NF_NOERR) THEN
+            PRINT*, "phyetat0: Lecture echouee pour <ALBLW>"
+            CALL abort
+         ENDIF
+         xmin = 1.0E+20
+         xmax = -1.0E+20
+         DO i = 1, klon
+            xmin = MIN(alblw(i,1),xmin)
+            xmax = MAX(alblw(i,1),xmax)
+         ENDDO
+         PRINT*,'Neige du sol <ALBLW>', xmin, xmax
+         DO nsrf = 2, nbsrf
+         DO i = 1, klon
+            alblw(i,nsrf) = alblw(i,1)
+         ENDDO
+         ENDDO
+      ENDIF
+
+cIM END alblw
+
+c
 c Lecture de evaporation:  
 c
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 467)
@@ -4,5 +4,5 @@
       SUBROUTINE phyredem (fichnom,dtime,radpas,
      .           rlat,rlon, pctsrf,tsol,tsoil,deltat,qsurf,qsol,snow,
-     .           albedo, evap, rain_fall, snow_fall,
+     .           albedo, alblw, evap, rain_fall, snow_fall,
      .           solsw, sollw,fder,
      .           radsol,frugs,agesno,
@@ -36,4 +36,7 @@
       REAL snow(klon,nbsrf)
       REAL albedo(klon,nbsrf)
+cIM BEG
+      REAL alblw(klon,nbsrf)
+cIM END
       REAL evap(klon,nbsrf)
       REAL rain_fall(klon)
@@ -360,4 +363,29 @@
 #endif
       ENDDO
+
+cIM BEG albedo LW
+        DO nsrf = 1, nbsrf
+        IF (nsrf.LE.99) THEN
+        WRITE(str2,'(i2.2)') nsrf
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid,"ALBLW"//str2,NF_DOUBLE,1,idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid,"ALBLW"//str2,NF_FLOAT,1,idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 23,
+     .                        "albedo LW de surface No."//str2)
+        ierr = NF_ENDDEF(nid)
+        ELSE
+        PRINT*, "Trop de sous-mailles"
+        CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,alblw(1,nsrf))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,alblw(1,nsrf))
+#endif
+      ENDDO
+cIM END albedo LW
 c
       DO nsrf = 1, nbsrf
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 467)
@@ -191,12 +191,158 @@
       PARAMETER(klevp1=klev+1)
 #include "raddim.h"
-      REAL*8 ZFSUP(KDLON,KFLEV+1)
-      REAL*8 ZFSDN(KDLON,KFLEV+1)
-      REAL*8 ZFSUP0(KDLON,KFLEV+1)
-      REAL*8 ZFSDN0(KDLON,KFLEV+1)
+cc      REAL*8 ZFSUP(KDLON,KFLEV+1)
+cc      REAL*8 ZFSDN(KDLON,KFLEV+1)
+cc      REAL*8 ZFSUP0(KDLON,KFLEV+1)
+cc      REAL*8 ZFSDN0(KDLON,KFLEV+1)
+c
+      REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2)
+      SAVE swdn0 , swdn, swup0, swup
 
 cccIM cf. FH
       real u850(klon),v850(klon),u200(klon),v200(klon)
       real u500(klon),v500(klon),phi500(klon),w500(klon)
+cIM
+      real prw(klon)
+
+cIM ISCCP - proprietes microphysiques des nuages convectifs
+      REAL convliq(klon,klev)  ! eau liquide nuageuse convective
+      REAL convfra(klon,klev)  ! fraction nuageuse convective
+
+      REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
+      REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
+      REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
+      REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
+
+      INTEGER kinv, linv
+
+cIM ISCCP simulator BEGIN
+      INTEGER igfi2D(iim,jjmp1)
+cv3.4
+      INTEGER debug, debugcol
+      INTEGER npoints
+      PARAMETER(npoints=klon) 
+      INTEGER sunlit(klon)
+
+      INTEGER ncol, seed(klon)
+
+cIM dans clesphys.h top_height, overlap
+c     PARAMETER(ncol=100)
+c     PARAMETER(ncol=625)
+      PARAMETER(ncol=10)
+      REAL tautab(0:255)
+      INTEGER invtau(-20:45000)
+      REAL emsfc_lw
+      PARAMETER(emsfc_lw=0.99)
+      REAL    ran0                      ! type for random number fuction
+
+      REAL pfull(klon,klev)
+      REAL phalf(klon,klev+1)
+      REAL cldtot(klon,klev)
+      REAL dtau_s(klon,klev)
+      REAL dtau_c(klon,klev)
+      REAL dem_s(klon,klev)
+      REAL dem_c(klon,klev)
+cPLUS : variables de haut en bas pour le simulateur ISCCP
+      REAL qv(klon,klev)
+      REAL cc(klon,klev)
+      REAL conv(klon,klev)
+      REAL dtau_sH2B(klon,klev)
+      REAL dtau_cH2B(klon,klev)
+      REAL at(klon,klev)
+      REAL dem_sH2B(klon,klev)
+      REAL dem_cH2B(klon,klev)
+
+c output from ISCCP
+      REAL fq_isccp(klon,7,7)
+      REAL totalcldarea(klon) 
+      REAL meanptop(klon)
+      REAL meantaucld(klon)
+      REAL boxtau(klon,ncol)
+      REAL boxptop(klon,ncol) 
+
+c grille 4d physique
+      INTEGER l, ni, nj, kmax, lmax, nrec
+      INTEGER ni1, ni2, nj1, nj2
+c     PARAMETER(kmax=7, lmax=7)
+      PARAMETER(kmax=8, lmax=8)
+      INTEGER kmaxm1, lmaxm1
+      PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
+c     INTEGER iimx7, jjmx7, jjmp1x7 
+c     PARAMETER(iimx7=iim*7, jjmx7=jjm*7, jjmp1x7=jjmp1*7)
+c     REAL fq4d(iim,jjmp1,7,7)
+c     REAL fq3d(iimx7, jjmp1x7)
+      INTEGER iimx8, jjmx8, jjmp1x8 
+      PARAMETER(iimx8=iim*8, jjmx8=jjm*8, jjmp1x8=jjmp1*8)
+      REAL fq4d(iim,jjmp1,8,8)
+      REAL fq3d(iimx8, jjmp1x8)
+cIM180603     SAVE fq3d
+
+c     REAL maxfq3d, minfq3d
+c
+      INTEGER iw, iwmax
+      REAL wmin, pas_w
+c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
+      PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
+      REAL o500(klon)
+      INTEGER nreg, nbreg
+      PARAMETER(nbreg=5)
+c     REAL histoW(iwmax,kmaxm1,lmaxm1)
+      REAL histoW(kmaxm1,lmaxm1,iwmax,nbreg)
+      REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbreg)
+cIM180603     
+c     SAVE histoW, nhistoW
+c     SAVE nhistoW
+      REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbreg) 
+      SAVE nhistoWt
+
+c     REAL histoWinv(kmaxm1,lmaxm1,iwmax)
+c     REAL nhistoW(kmaxm1,lmaxm1,iwmax)
+      INTEGER linv
+c     LOGICAL pct_ocean(klon,nbreg)
+      INTEGER pct_ocean(klon,nbreg)
+      REAL rlonPOS(klon) 
+c     CHARACTER*4 pdirect
+ 
+c sorties ISCCP
+
+      logical ok_isccp
+      real ecrit_isccp
+      integer nid_isccp
+      save ok_isccp, ecrit_isccp, nid_isccp        
+
+#define histISCCP
+#undef histISCCP
+#ifdef histISCCP
+c     data ok_isccp,ecrit_isccp/.true.,0.125/      
+c     data ok_isccp,ecrit_isccp/.true.,1./      
+      data ok_isccp/.true./      
+#else
+      data ok_isccp/.false./
+#endif
+
+      REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
+c     DATA zx_tau/0.1, 1.3, 3.6, 9.4, 23., 60./
+c     DATA zx_pc/50., 180., 310., 440., 560., 680., 800., 1015./
+c     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
+cOK     DATA zx_tau/0.0, 0.1, 1.3, 3.6, 9.4, 23., 60./
+cOK     DATA zx_pc/800., 680., 560., 440., 310., 180., 50./
+
+c tester l'alure
+      DATA zx_tau/1., 2., 3., 4., 5., 6., 7./
+c     DATA zx_pc/1., 2., 3., 4., 5., 6., 7./
+      DATA zx_pc/7., 6., 5., 4., 3., 2., 1./
+
+      INTEGER komega, nhoriRD 
+
+c statistiques regime dynamique END
+
+c     REAL del_lon(iim), del_lat(jjmp1)
+      REAL del_lon, del_lat
+c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
+      REAL zx_lonx8(iimx8), zx_latx8(jjmp1x8)
+c     INTEGER nhorix7
+      INTEGER nhorix8
+
+cIM ISCCP simulator END
 
       logical ok_hf
@@ -497,5 +643,4 @@
       SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
 cccIM
-      SAVE  ZFSUP,ZFSDN,ZFSUP0,ZFSDN0 
 
       INTEGER itaprad
@@ -753,5 +898,5 @@
          CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
      .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow,
-     .       falbe, fevap, rain_fall,snow_fall,solsw, sollwdown,
+     .       falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
      .       dlw,radsol,frugs,agesno,clesphy0,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
@@ -880,4 +1025,9 @@
 c   Initialisation des sorties
 c=============================================================
+
+#ifdef histISCCP
+#include "ini_histISCCP.h"
+#endif
+
 #ifdef histhf
 #include "ini_histhf.h"
@@ -958,5 +1108,4 @@
       ENDIF
 C
-      IF (if_ebil.ge.1) THEN 
         DO i = 1, klon
           ztsol(i) = 0.
@@ -967,4 +1116,6 @@
           ENDDO
         ENDDO
+C
+      IF (if_ebil.ge.1) THEN 
         ztit='after dynamic'
         CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
@@ -1072,6 +1223,6 @@
       DO nsrf = 1, nbsrf
       DO i = 1, klon
-         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
-cccc        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
+c         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
+        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
       ENDDO
       ENDDO
@@ -1091,5 +1242,14 @@
         rmu0 = -999.999
       ENDIF
-C
+cIM BEG
+      DO i=1, klon
+       sunlit(i)=1 
+       IF(rmu0(i).EQ.0.) sunlit(i)=0
+c      IF(rmu0(i).EQ.0.) THEN 
+c       sunlit(i)=0 
+c       PRINT*,' il fait nuit ',i,rlat(i),rlon(i)
+c      ENDIF
+      ENDDO
+cIM END
 C     Calcul de l'abedo moyen par maille
       albsol(:)=0.
@@ -1103,10 +1263,16 @@
 C
 C     Repartition sous maille des flux LW et SW
-      DO nsrf = 1, nbsrf
-      DO i = 1, klon
-        fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4
-        fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
-      ENDDO
-      ENDDO
+C Modif OM+PASB+JLD
+C Repartition du longwave par sous-surface linearisee
+Cn
+       DO nsrf = 1, nbsrf
+       DO i = 1, klon
+c$$$        fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4
+c$$$        fsollw(i,nsrf) = sollw(i)
+         fsollw(i,nsrf) = sollw(i)
+     $      + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i,nsrf))
+         fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
+       ENDDO
+       ENDDO
 
       fder = dlw
@@ -1116,5 +1282,5 @@
      e            julien, rmu0, 
      e            ok_veget, ocean, npas, nexca, ftsol,
-     $            soil_model,ftsoil, qsol, 
+     $            soil_model,cdmmax, cdhmax, ftsoil, qsol, 
      $            paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
      $            fluxlat,
@@ -1616,4 +1782,574 @@
       enddo
 
+cIM ISCCP simulator BEGIN
+      IF (ok_isccp) THEN
+cIM calcul tau. emi nuages convectifs
+      convfra(:,:)=rnebcon(:,:)
+      convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
+c     CALL newmicro (paprs, pplay,ok_newmicro,
+c    .            t_seri, cldliq, cldfra, cldtau, cldemi,
+c    .            cldh, cldl, cldm, cldt, cldq)
+      CALL newmicro (paprs, pplay,ok_newmicro,
+     .            t_seri, convliq, convfra, dtau_c, dem_c,
+     .            cldh_c, cldl_c, cldm_c, cldt_c, cldq_c)
+
+cIM calcul tau. emi nuages startiformes
+      CALL newmicro (paprs, pplay,ok_newmicro,
+     .            t_seri, cldliq, cldfra, dtau_s, dem_s,
+     .            cldh_s, cldl_s, cldm_s, cldt_s, cldq_s)
+cIM calcul diagramme (PC, tau) cf. ISCCP D
+c     seed=50
+c     seed=ran0(klon)
+cT1O3     
+c     top_height=1
+cT3O3
+c     top_height=3
+c     overlap=3
+cIM cf GCM     
+      cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
+
+cIM inversion des niveaux de pression ==> de haut en bas
+      DO k=1,klev
+       kinv=klev-k+1
+       DO i=1,klon
+        pfull(i,k)=pplay(i,kinv)
+c on met toutes les variables de Haut 2 Bas
+        qv(i,k)=q_seri(i,kinv)
+        cc(i,k)=cldtot(i,kinv)
+        conv(i,k)=rnebcon(i,kinv)
+        dtau_sH2B(i,k)=dtau_s(i,kinv)
+        dtau_cH2B(i,k)=dtau_c(i,kinv)
+        at(i,k)=t_seri(i,kinv)
+        dem_sH2B(i,k)=dem_s(i,kinv)
+        dem_cH2B(i,k)=dem_c(i,kinv)
+
+       ENDDO 
+      ENDDO 
+
+      DO k=1,klev+1
+       kinv=klev-k+2
+       DO i=1,klon
+        phalf(i,k)=paprs(i,kinv)
+       ENDDO 
+      ENDDO 
+
+c     open(99,file='tautab.bin',access='sequential',
+c    $     form='unformatted',status='old')
+c     read(99) tautab
+
+cIM210503
+      IF (debut) THEN
+      open(99,file='tautab.formatted', FORM='FORMATTED')
+      read(99,'(f30.20)') tautab
+      close(99)
+c
+      open(99,file='invtau.formatted',form='FORMATTED')
+      read(99,'(i10)') invtau
+      close(99)
+c
+       nsrf=3
+       DO nreg=1, nbreg
+       DO i=1, klon
+
+c       IF (debut) THEN
+         IF(rlon(i).LT.0.) THEN
+           rlonPOS(i)=rlon(i)+360.
+         ELSE
+           rlonPOS(i)=rlon(i)  
+         ENDIF
+c       ENDIF
+
+c       pct_ocean(i,nreg)=.FALSE.
+        pct_ocean(i,nreg)=0
+
+c      DO nsrf = 1, nbsrf 
+
+c test si c'est 1 point d'ocean
+        IF(pctsrf(i,nsrf).EQ.1.) THEN
+
+         IF(nreg.EQ.1) THEN
+
+c TROP
+          IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
+c          pct_ocean(i,nreg)=.TRUE.
+           pct_ocean(i,nreg)=1
+          ENDIF
+
+c PACIFIQUE NORD
+          ELSEIF(nreg.EQ.2) THEN
+           IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
+            IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN 
+c            pct_ocean(i,nreg)=.TRUE.
+             pct_ocean(i,nreg)=1
+            ENDIF
+           ENDIF
+c CALIFORNIE ST-CU
+         ELSEIF(nreg.EQ.3) THEN
+          IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN
+           IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
+c           pct_ocean(i,nreg)=.TRUE.
+            pct_ocean(i,nreg)=1
+           ENDIF
+          ENDIF
+c HAWAI
+        ELSEIF(nreg.EQ.4) THEN 
+         IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
+          IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
+c          pct_ocean(i,nreg)=.TRUE.
+           pct_ocean(i,nreg)=1
+          ENDIF
+         ENDIF
+c WARM POOL
+        ELSEIF(nreg.EQ.5) THEN 
+         IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
+          IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
+c          pct_ocean(i,nreg)=.TRUE.
+           pct_ocean(i,nreg)=1
+          ENDIF
+         ENDIF
+        ENDIF !nbreg
+c TROP
+c        IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
+c         pct_ocean(i)=.TRUE.
+c         WRITE(*,*) 'pct_ocean =',i, rlon(i), rlat(i)
+c          ENDIF !lon
+c         ENDIF !lat
+
+        ENDIF !pctsrf
+c      ENDDO 
+       ENDDO !klon
+       ENDDO !nbreg
+ 
+cIM somme de toutes les nhistoW BEG
+      DO nreg = 1, nbreg
+      DO k = 1, kmaxm1
+      DO l = 1, lmaxm1
+      DO iw = 1, iwmax
+       nhistoWt(k,l,iw,nreg)=0.
+      ENDDO
+      ENDDO
+      ENDDO
+      ENDDO
+cIM somme de toutes les nhistoW END
+      ENDIF
+
+
+c     CALL ISCCP_CLOUD_TYPES(nlev,ncol,seed,pfull,phalf,qv,
+c    &     cc,conv,dtau_s,dtau_c,top_height,overlap,
+c    &     tautab,invtau,skt,emsfc_lw,at,dem_s,dem_c,fq_isccp,
+c    &     totalcldarea,meanptop,meantaucld,boxtau,boxptop)
+
+c     DO i=1, klon
+c     i=1
+c1011  CONTINUE
+c
+cIM on verifie les donnees de INPUT en dehors du simulateur ISCCP
+cIM 1D non-vectorise (!) pour qu'on gagne du temps ...
+cIM
+c BEGIN find unpermittable data.....
+!     ---------------------------------------------------!
+!     find unpermittable data.....
+!
+c     do 13 k=1,klev
+c ca prend trop de temps ??
+c     cldtot(:,:) = min(max(cldtot(:,:),0.),1.)
+c     rnebcon(:,:) = min(max(rnebcon(:,:),0.),1.)
+c     dtau_s(:,:) = max(dtau_s(:,:),0.)
+c     dem_s(:,:) = min(max(dem_s(:,:),0.),1.)
+c     dtau_c(:,:) = max(dtau_c(:,:),0.)
+c     dem_c(:,:) = min(max(dem_c(:,:),0.),1.)
+c ca prend trop de temps ??
+
+c           if (cldtot(i,k) .lt. 0.) then
+c               print *, ' error = cloud fraction less than zero'
+c               STOP
+c           end if
+c           if (cldtot(i,k) .gt. 1.) then
+c               print *, ' error = cloud fraction greater than 1'
+c               STOP
+c           end if
+c           if (rnebcon(i,k) .lt. 0.) then
+c               print *,
+c    &           ' error = convective cloud fraction less than zero'
+c               STOP
+c           end if
+c           if (rnebcon(i,k) .gt. 1.) then
+c               print *,
+c    &           ' error = convective cloud fraction greater than 1'
+c               STOP
+c           end if
+
+c           if (dtau_s(i,k) .lt. 0.) then
+c               print *,
+c    &           ' error = stratiform cloud opt. depth less than zero'
+c               STOP
+c           end if
+c           if (dem_s(i,k) .lt. 0.) then
+c               print *,
+c    &           ' error = stratiform cloud emissivity less than zero'
+c               STOP
+c           end if
+c           if (dem_s(i,k) .gt. 1.) then
+c               print *,
+c    &           ' error = stratiform cloud emissivity greater than 1'
+c               STOP
+c           end if
+
+c           if (dtau_c(i,k) .lt. 0.) then
+c               print *,
+c    &           ' error = convective cloud opt. depth less than zero'
+c               STOP
+c           end if
+c           if (dem_c(i,k) .lt. 0.) then
+c               print *,
+c    &           ' error = convective cloud emissivity less than zero'
+c               STOP
+c           end if
+c           if (dem_c(i,k) .gt. 1.) then
+c               print *,
+c    &           ' error = convective cloud emissivity greater than 1'
+c               STOP
+c           end if
+c13    continue
+
+!     ---------------------------------------------------!
+c
+c END   find unpermittable data.....
+cv2.2.1.1     DO i=1, klon
+c     i=1
+c     seed=i
+c
+cv3.4
+      if (debut) then
+        DO i=1, klon
+          seed(i)=i+100
+c         seed(i)=i+50
+        ENDDO
+      endif
+c     seed=aint(ran0(klon))
+c     CALL ISCCP_CLOUD_TYPES(klev,ncol,seed,pfull(i,:),phalf(i,:)
+cv2.2.1.1
+c     CALL ISCCP_CLOUD_TYPES(klev,ncol,seed(i),pfull(i,:),phalf(i,:)
+c    &     ,q_seri(i,:),
+c    &     cldtot(i,:),rnebcon(i,:),dtau_s(i,:),dtau_c(i,:),
+c    &     top_height,overlap,
+c    &     tautab,invtau,ztsol,emsfc_lw,t_seri(i,:),dem_s(i,:),
+c    &     dem_c(i,:),
+c    &     fq_isccp(i,:,:),
+c    &     totalcldarea(i),meanptop(i),meantaucld(i),
+c    &     boxtau(i,:),boxptop(i,:))
+cv2.2.1.1
+cv3.4
+      debug=0
+      debugcol=0
+cIM260503
+c o500 ==> distribution nuage ftion du regime dynamique 
+      DO i=1, klon
+       o500(i)=omega(i,8)*864.
+c      PRINT*,'pphi8 ',pphi(i,8),'zphi8,11,12',zphi(i,8),
+c    & zphi(i,11),zphi(i,12)
+      ENDDO
+
+c axe vertical pour les differents niveaux des histogrammes
+c     DO iw=1, iwmax
+c       zx_o500(iw)=wmin+(iw-1./2.)*pas_w
+c     ENDDO 
+c     PRINT*,' phys AVANT seed(3361)=',seed(3361) 
+      CALL ISCCP_CLOUD_TYPES(
+     &     debug,
+     &     debugcol,
+     &     klon,
+     &     sunlit,
+     &     klev,
+     &     ncol,
+     &     seed,
+     &     pfull,
+     &     phalf,
+c var de bas en haut ==> PB !
+c    &     q_seri,
+c    &     cldtot,
+c    &     rnebcon,
+c    &     dtau_s,
+c    &     dtau_c,
+c var de Haut en Bas BEG
+     &     qv, cc, conv, dtau_sH2B, dtau_cH2B,
+c var de Haut en Bas END
+     &     top_height,
+     &     overlap,
+     &     tautab,
+     &     invtau,
+     &     ztsol,
+     &     emsfc_lw,
+c var de bas en haut ==> PB !
+c    &     t_seri,
+c    &     dem_s,
+c    &     dem_c,
+c var de Haut en Bas BEG
+     &     at, dem_sH2B, dem_cH2B,
+cIM260503
+c    &     o500, pct_ocean,
+c var de Haut en Bas END
+     &     fq_isccp,
+     &     totalcldarea,
+     &     meanptop,
+     &     meantaucld,
+     &     boxtau,
+     &     boxptop)
+c    &     boxptop,
+cIM 260503
+c    &     histoW,
+c    &     nhistoW    
+c    &)
+
+cIM 200603
+c     PRINT*,'physiq fq_isccp(6,1,1)',fq_isccp(6,1,1)
+       
+cIM 200603
+cIM somme de toutes les nhistoW BEG
+c     DO k = 1, kmaxm1
+c     DO l = 1, lmaxm1
+c     DO iw = 1, iwmax
+c     nhistoWt(k,l,iw)=nhistoWt(k,l,iw)+nhistoW(k,l,iw)
+ccc      IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
+c      IF(nhistoWt(k,l,iw).NE.0.) THEN 
+c       PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
+c      ENDIF
+c     ENDDO
+c     ENDDO
+c     ENDDO
+cIM somme de toutes les nhistoW END
+c     PRINT*,' phys APRES seed(3361)=',seed(3361) 
+cv3.4
+c     i=i+1
+c     IF(i.LE.klon) THEN
+c      GOTO 1011
+c     ENDIF
+cv2.2.1.1     ENDDO 
+
+c passage de la grille (klon,7,7) a (iim,jjmp1,7,7)
+c     minfq3d=100.
+c     maxfq3d=0.
+cIM calcul des correspondances entre la grille physique et
+cIM la grille dynamique
+c     DO i=1, klon
+c       grid_phys(i)=i
+c       PRINT*,'i, grid_phys',i,grid_phys(i)
+c     ENDDO
+c     CALL gr_fi_dyn(1,klon,iimp1,jjmp1,grid_phys,grid_dyn)
+c     DO j=1, jjmp1 
+c       DO i=1, iimp1
+c        PRINT*,'i,j grid_dyn ',i,j,grid_dyn(i,j)
+c       ENDDO
+c     ENDDO
+c
+      DO l=1, lmax
+       DO k=1, kmax
+cIM grille physique ==> grille ecriture 2D (iim,jjmp1)
+c
+        DO i=1, iim
+          fq4d(i,1,k,l)=fq_isccp(1,k,l) 
+cc         PRINT*,'first j=1 i =',i
+        ENDDO
+        DO j=2, jjm
+          DO i=1, iim
+cERROR ??         ig=i+iim*(j-1)
+          ig=i+1+(j-2)*iim
+cc         PRINT*,'i =',i,'j =',j,'ig=',ig
+          fq4d(i,j,k,l)=fq_isccp(ig,k,l)              
+         ENDDO
+        ENDDO
+        DO i=1, iim
+          fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l) 
+cc         PRINT*,'last jjmp1 i =',i
+        ENDDO
+        IF(debut) THEN
+        DO j=1, jjmp1
+          DO i=1, iim
+            IF(j.GE.2.AND.j.LE.jjm) THEN
+              igfi2D(i,j)=i+1+(j-2)*iim
+c             PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
+            ELSEIF(j.EQ.1) THEN
+              igfi2D(i,j)=1
+c             PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
+            ELSEIF(j.EQ.jjmp1) THEN
+              igfi2D(i,j)=klon
+c             PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
+            ENDIF
+          ENDDO
+        ENDDO
+        ENDIF
+c       STOP
+c
+c       CALL gr_fi_ecrit(1,klon,iim,jjmp1,fq_isccp(:,k,l),
+c    $       fq4d(:,:,k,l))
+       ENDDO 
+      ENDDO 
+      DO l=1, lmax
+       fq4d(:,:,8,l)=-1.e+10
+       fq4d(:,:,l,8)=-1.e+10
+      ENDDO
+      DO l=1, lmax
+       DO k=1, kmax  
+        DO j=1, jjmp1
+         DO i=1, iim
+
+c inversion TAU ?!
+c         ni=(i-1)*lmax+l
+c         nj=(j-1)*kmax+kmax-k+1
+c
+c210503 inversion en PC ==> pas besoin !!!
+c         ni=(i-1)*lmax+lmax-l+1
+c         nj=(j-1)*kmax+k
+c
+c210503
+          ni=(i-1)*lmax+l
+          nj=(j-1)*kmax+k
+ 
+c210503         if(k.EQ.8) then 
+c          fq4d(i,j,8,l)=-1.e+10
+c         endif
+
+c210503         if(l.EQ.8) then 
+c          fq4d(i,j,k,8)=-1.e+10
+c         endif
+
+          fq3d(ni,nj)=fq4d(i,j,k,l)
+
+c         if(fq3d(ni,nj).lt.0.) then
+c          print*,' fq3d LT ZERO ',ni,nj,fq3d(ni,nj)
+c         endif
+c         if(fq3d(ni,nj).gt.100.) then
+c          print*,' fq3d GT 100 ',ni,nj,fq3d(ni,nj)
+c         endif
+c max & min fq3d
+c         if(fq3d(ni,nj).gt.maxfq3d) maxfq3d=fq3d(ni,nj)
+c         if(fq3d(ni,nj).lt.minfq3d) minfq3d=fq3d(ni,nj)
+         
+         ENDDO
+        ENDDO
+c       fq4d(:,:,8,l)=-1.e+10
+c       fq4d(:,:,k,8)=-1.e+10
+c       k=k+1
+c       if(k.LE.kmax) then
+c        goto 1022
+c       endif
+       ENDDO
+c      l=l+1
+c      if(l.LE.lmax) then
+c       goto 1021
+c      endif
+      ENDDO
+
+c     print*,' minfq3d=',minfq3d,' maxfq3d=',maxfq3d
+c
+c calculs statistiques distribution nuage ftion du regime dynamique 
+c     DO i=1, klon
+c!      o500(i)=omega(i,9)*864.
+c!      PRINT*,' o500=',o500(i),' pphi(9)=',pphi(i,9)
+c       o500(i)=omega(i,8)*864.
+cc      PRINT*,' pphi(8)',pphi(i,8),'pphi(11)',pphi(i,11),
+cc    .'pphi(12)',pphi(i,12)
+cc      PRINT*,' zphi8,11,12=',zphi(i,8),zphi(i,11),zphi(i,12)
+cc     PRINT*,' o500',o500(i),' w500',w500(i)
+c     ENDDO
+
+c axe vertical pour les differents niveaux des histogrammes
+c     DO iw=1, iwmax
+c       zx_o500(iw)=wmin+(iw-1./2.)*pas_w
+c     ENDDO 
+
+
+c Ce calcul doit etre fait a partir de valeurs mensuelles ??
+cc     CALL histo_o500_pctau(o500,fq4d,histoW)
+cc     CALL histo_o500_pctau(paire,pctsrf,o500,fq4d,histoW)
+cc     CALL histo_o500_pctau(pct_ocean,rlat,o500,fq4d,histoW)
+ccOK ???     CALL histo_o500_pctau(pct_ocean,o500,fq4d,histoW)
+c     CALL histo_o500_pctau(klon,pct_ocean,o500,fq4d,histoW,nhistoW)
+c     CALL histo_o500_pctau(klon,pct_ocean,o500,fq_isccp,
+      CALL histo_o500_pctau(nbreg,pct_ocean,o500,fq_isccp,
+     &histoW,nhistoW)
+c
+cIM somme de toutes les nhistoW BEG
+      DO nreg=1, nbreg
+      DO k = 1, kmaxm1
+      DO l = 1, lmaxm1
+      DO iw = 1, iwmax
+       nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+
+     & nhistoW(k,l,iw,nreg)
+ccc      IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
+c      IF(nhistoWt(k,l,iw).NE.0.) THEN 
+c       PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
+c      ENDIF
+      ENDDO
+      ENDDO
+      ENDDO
+      ENDDO
+cIM somme de toutes les nhistoW END
+c
+c     IF(lafin) THEN   
+c     DO nreg=1, nbreg
+c     DO iw=1, iwmax
+c     DO l=1,lmaxm1
+c     DO k=1,kmaxm1
+c      IF(histoW(k,l,iw,nreg).NE.0.) then
+c        PRINT*,'physiq H nH',k,l,iw,
+c    &       histoW(k,l,iw,nreg),
+c    &       nhistoW(k,l,iw,nreg),nhistoWt(k,l,iw,nreg)
+c      ENDIF
+c     ENDDO
+c     ENDDO
+c     ENDDO
+c     ENDDO
+cIM verif fq_isccp, fq4d, fq3d
+c     DO l=1, lmaxm1
+c       DO k=1,kmaxm1
+c     i=74
+c     j=36
+c     DO j=1, jjmp1
+c      DO i=1, iim
+c       DO l=1, lmaxm1
+c         WRITE(*,'(a,3i4,7f10.4)')
+c    &    'fq_isccp,j,i,l=',j,i,l,
+c    &    (fq_isccp(igfi2D(i,j),k,l),k=1,kmaxm1)
+c         WRITE(*,'(a,3i4,7f10.4)')
+c    &    'fq4d,j,i,l=',j,i,l,(fq4d(i,j,k,l),k=1,kmaxm1)
+c       ENDDO
+c      ENDDO
+c     ENDDO
+c     ni1=(i-1)*8+1
+c     ni2=i*8
+c     nj1=(j-1)*8+1
+c     nj2=j*8
+c     DO ni=ni1,ni2
+c     WRITE(*,'(a,2i4,7f10.4)')
+c    &     'fq3d, ni,nj=',ni,nj,
+c    &      (fq3d(ni,nj),nj=nj1,nj2)
+c     ENDDO
+c     ENDIF
+
+c     DO iw=1, iwmax
+c      DO l=1,lmaxm1
+c       DO k=1,kmaxm1
+c        PRINT*,' iw,l,k,nhistoW=',iw,l,k,nhistoW(k,l,iw)
+c       ENDDO
+c      ENDDO
+c     ENDDO
+
+c       DO iw=1, iwmax
+c        DO l=1, lmaxm1
+c         linv=lmaxm1-l+1
+c         DO k=1, kmaxm1
+c         histoWinv(k,l,iw)=histoW(iw,k,l)
+c       ENDDO
+c      ENDDO
+c     ENDDO
+c
+c pb syncronisation ?? : 48 * 30 * 7 (jour1) + 48* 29 * 7 (jour suivant)
+c
+
+
+      ENDIF !ok_isccp
+cIM ISCCP simulator END
+
 c   On prend la somme des fractions nuageuses et des contenus en eau
       cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
@@ -1717,5 +2453,5 @@
 cccIMs             topsw0,toplw0,solsw0,sollw0)
      s             topsw0,toplw0,solsw0,sollw0,
-     s             ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
+     s             swdn0, swdn, swup0, swup     )
       itaprad = 0
       ENDIF
@@ -1968,4 +2704,16 @@
 cIM cf. FH     slp(:) = paprs(:,1)*exp(pphis(:)/(289.*t_seri(:,1)))
       slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
+c     PRINT*,' physiq slp ',slp(2185),paprs(2185,1),pphis(2185),
+c    .       RD,t_seri(2185,1)
+c
+ccc prw = eau precipitable
+      DO i = 1, klon
+       prw(i) = 0.
+      DO k = 1, klev
+       prw(i) = prw(i) +
+     .          q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
+      ENDDO
+c      PRINT*,' i ',i,' prw',prw(i)
+      ENDDO
 c
 
@@ -1973,4 +2721,8 @@
 c   Ecriture des sorties
 c=============================================================
+
+#ifdef histISCCP
+#include "write_histISCCP.h"
+#endif
 
 #ifdef histhf
@@ -2024,5 +2776,5 @@
          CALL phyredem ("restartphy.nc",dtime,radpas,
      .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
-     .      fsnow, falbe, fevap, rain_fall, snow_fall,
+     .      fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
      .      solsw, sollwdown,dlw,
      .      radsol,frugs,agesno,
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 467)
@@ -7,5 +7,5 @@
      .                  sollwdown,
      .                  topsw0,toplw0,solsw0,sollw0,
-     .                  ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
+     .                  swdn0, swdn, swup0, swup    )
       IMPLICIT none
 c======================================================================
@@ -98,4 +98,5 @@
       REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
       REAL*8 zznormcp
+      REAL swdn(klon,2),swdn0(klon,2),swup(klon,2),swup0(klon,2)
 c
 c-------------------------------------------
@@ -234,4 +235,12 @@
          sollw0(iof+i) = zsollw0(i)
          albpla(iof+i) = zalbpla(i)
+         swdn0 ( iof+i,1)   = ZFSDN0 ( i,1 )
+         swdn0 ( iof+i,2)   = ZFSDN0 ( i,kflev + 1 )
+         swdn  ( iof+i,1)   = ZFSDN  ( i,1 )
+         swdn  ( iof+i,2)   = ZFSDN  ( i,kflev + 1 )
+         swup0 ( iof+i,1)   = ZFSUP0 ( i,1 )
+         swup0 ( iof+i,2)   = ZFSUP0 ( i,kflev + 1 )
+         swup  ( iof+i,1)   = ZFSUP  ( i,1 )
+         swup  ( iof+i,2)   = ZFSUP  ( i,kflev + 1 )
       ENDDO
       DO k = 1, kflev
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histday.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histday.h	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histday.h	(revision 467)
@@ -160,24 +160,63 @@
 c    .              'Cloud liquid water path','-')
 c
-      zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, klevp1)
+c     zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, klevp1)
+      zx_tmp_fi2d(1 : klon) = swup( 1 : klon, 2 )
+
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_day, "SWupTOA",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, 1)
+c      zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swup( 1 : klon, 1 )
+
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_day, "SWupSFC",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, klevp1)
+c      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, klevp1)
+      zx_tmp_fi2d(1 : klon) = swdn( 1 : klon, 2 )
+
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_day, "SWdnTOA",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, 1)
+c      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swdn( 1 : klon, 1 )
+
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_day, "SWdnSFC",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 
+c      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, klevp1)
+      zx_tmp_fi2d(1 : klon) = swup0( 1 : klon, 2 )
+
+      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
+      CALL histwrite(nid_day, "SWupTOAclr",itau_w,zx_tmp_2d,
+     .                               iim*jjmp1,ndex2d)
+
+c      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swup0( 1 : klon, 1 )
+
+      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
+      CALL histwrite(nid_day, "SWupSFCclr",itau_w,zx_tmp_2d,
+     .                               iim*jjmp1,ndex2d)
+ 
+c      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, klevp1)
+      zx_tmp_fi2d(1 : klon) = swdn0( 1 : klon, 2 )
+
+      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
+      CALL histwrite(nid_day, "SWdnTOAclr",itau_w,zx_tmp_2d,
+     .                               iim*jjmp1,ndex2d)
+ 
+c      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swdn0( 1 : klon, 1 )
+
+      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
+      CALL histwrite(nid_day, "SWdnSFCclr",itau_w,zx_tmp_2d,
+     .                               iim*jjmp1,ndex2d)
+cIM
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, prw,zx_tmp_2d)
+      CALL histwrite(nid_day,"prw",itau_w,zx_tmp_2d,
+     .               iim*jjmp1,ndex2d)
 
 c   Ecriture de champs dynamiques sur des niveaux de pression
@@ -211,23 +250,4 @@
      .               iim*jjmp1,ndex2d)
 
-      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, klevp1)
-      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_day, "SWupTOAclr",itau_w,zx_tmp_2d,
-     .                               iim*jjmp1,ndex2d)
-
-      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, 1)
-      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_day, "SWupSFCclr",itau_w,zx_tmp_2d,
-     .                               iim*jjmp1,ndex2d)
- 
-      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, klevp1)
-      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_day, "SWdnTOAclr",itau_w,zx_tmp_2d,
-     .                               iim*jjmp1,ndex2d)
- 
-      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, 1)
-      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_day, "SWdnSFCclr",itau_w,zx_tmp_2d,
-     .                               iim*jjmp1,ndex2d)
 
 c
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h	(revision 466)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h	(revision 467)
@@ -459,41 +459,52 @@
 c
 cccIM
-      zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, klevp1)
+c      zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, klevp1)
+      zx_tmp_fi2d(1 : klon) = swup ( 1 : klon, 2 )
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWupTOA",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, 1)
+c      zx_tmp_fi2d(1 : klon) = ZFSUP( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swup ( 1 : klon, 1 )
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWupSFC",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, klevp1)
+c      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, klevp1)
+      zx_tmp_fi2d(1 : klon) = swdn ( 1 : klon, 2 )
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWdnTOA",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, 1)
+c      zx_tmp_fi2d(1 : klon) = ZFSDN( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swdn ( 1 : klon, 1 )
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWdnSFC",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, prw,zx_tmp_2d)
+      CALL histwrite(nid_mth,"prw",itau_w,zx_tmp_2d,
+     .               iim*jjmp1,ndex2d)
+
 cccIM clear sky
-      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, klevp1)
+c      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, klevp1)
+      zx_tmp_fi2d(1 : klon) = swup0 ( 1 : klon, 2 )
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWupTOAclr",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, 1)
+c      zx_tmp_fi2d(1 : klon) = ZFSUP0( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swup0 ( 1 : klon, 1 )
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWupSFCclr",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, klevp1)
+c      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, klevp1)
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWdnTOAclr",itau_w,zx_tmp_2d,
      .                               iim*jjmp1,ndex2d)
 c
-      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, 1)
+c      zx_tmp_fi2d(1 : klon) = ZFSDN0( 1 : klon, 1)
+      zx_tmp_fi2d(1 : klon) = swdn0 ( 1 : klon, 1 )
       CALL gr_fi_ecrit(1,klon,iim,jjmp1,zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth, "SWdnSFCclr",itau_w,zx_tmp_2d,
