Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/albedo.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/albedo.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/albedo.F	(revision 486)
@@ -182,32 +182,2 @@
       END
 c========================================================================
-      SUBROUTINE albsno(veget, agesno, alb_neig)
-      IMPLICIT none
-c
-#include "dimensions.h"
-#include "dimphy.h"
-      INTEGER nvm
-      PARAMETER (nvm=8)
-      REAL veget(klon,nvm)
-      REAL alb_neig(klon)
-      REAL agesno(klon)
-c
-      INTEGER i, nv
-c
-      REAL init(nvm), decay(nvm), as
-      SAVE init, decay
-      DATA init /0.55, 0.14, 0.18, 0.29, 0.15, 0.15, 0.14, 0./
-      DATA decay/0.30, 0.67, 0.63, 0.45, 0.40, 0.14, 0.06, 1./
-c
-      DO i = 1, klon
-         alb_neig(i) = 0.0
-      ENDDO
-      DO nv = 1, nvm
-         DO i = 1, klon
-            as = init(nv)+decay(nv)*EXP(-agesno(i)/5.)
-            alb_neig(i) = alb_neig(i) + veget(i,nv)*as
-         ENDDO
-      ENDDO
-c
-      RETURN
-      END
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.h	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.h	(revision 486)
@@ -11,4 +11,8 @@
 cIM seuils cdrm, cdrh
        REAL cdmmax, cdhmax
+cIM param. stabilite s/ terres et en dehors
+       REAL ksta, ksta_ter
+cIM ok_kzmin : clef calcul Kzmin dans la CL de surface cf FH
+       LOGICAL ok_kzmin
 
        COMMON/clesphys/cycle_diurne, soil_model, new_oliq,
@@ -16,3 +20,4 @@
      ,     , co2_ppm, solaire, RCO2, RCH4, RN2O, RCFC11, RCFC12
      ,     , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt
-     ,     , top_height, overlap, cdmmax, cdhmax
+     ,     , top_height, overlap, cdmmax, cdhmax, ksta, ksta_ter
+     ,     , ok_kzmin
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.inc
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.inc	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clesphys.inc	(revision 486)
@@ -9,4 +9,6 @@
        INTEGER :: top_height, overlap
        REAL :: cdmmax, cdhmax
+       REAL :: ksta, ksta_ter
+       LOGICAL :: ok_kzmin
 
        COMMON/clesphys/cycle_diurne, soil_model, new_oliq, &
@@ -14,3 +16,4 @@
      &     , co2_ppm, solaire, RCO2, RCH4, RN2O, RCFC11, RCFC12 &
      &     , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt &
-     &     , top_height, overlap, cdmmax, cdhmax
+     &     , top_height, overlap, cdmmax, cdhmax, ksta, ksta_ter &
+     &     , ok_kzmin
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 486)
@@ -7,5 +7,6 @@
      .                  jour, rmu0, co2_ppm,
      .                  ok_veget, ocean, npas, nexca, ts,
-     .                  soil_model,cdmmax, cdhmax, ftsoil,qsol,
+     .                  soil_model,cdmmax, cdhmax,
+     .                  ksta, ksta_ter, ok_kzmin, ftsoil,qsol,
      .                  paprs,pplay,radsol,snow,qsurf,evap,albe,alblw,
      .                  fluxlat,
@@ -143,4 +144,8 @@
 cIM ajout seuils cdrm, cdrh
       REAL cdmmax, cdhmax
+cIM: 261103
+      REAL ksta, ksta_ter
+      LOGICAL ok_kzmin
+cIM: 261103
       REAL ftsoil(klon,nsoilmx,nbsrf)
       REAL ytsoil(klon,nsoilmx)
@@ -472,4 +477,7 @@
 c calculer Cdrag et les coefficients d'echange
       CALL coefkz(nsrf, knon, ypaprs, ypplay,
+cIM 261103
+     .     ksta, ksta_ter,
+cIM 261103
      .            yts, yrugos, yu, yv, yt, yq,
      .            yqsurf, 
@@ -494,4 +502,27 @@
       endif
 
+c
+cIM: 261103
+      if (ok_kzmin) THEN
+cIM cf FH: 201103 BEG
+c   Calcul d'une diffusion minimale pour les conditions tres stables.
+      call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm
+     .   ,ycoefm0,ycoefh0)
+c      call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ         ')
+c      call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN      ')
+ 
+       if ( 1.eq.1 ) then
+       DO k = 1, klev
+       DO i = 1, knon
+          ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
+          ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
+       ENDDO
+       ENDDO
+       endif
+cIM cf FH: 201103 END
+      endif !ok_kzmin
+cIM: 261103
+
+c
 c calculer la diffusion des vitesses "u" et "v"
       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
@@ -503,11 +534,4 @@
       ytaux = y_flux_u(:,1)
       ytauy = y_flux_v(:,1)
-
-      if (nsrf.eq.is_oce) then
-         do j=1,knon
-            ycoefm(j,1)=1.e-3
-c            ycoefh(j,1)=0.8*ycoefm(j,1)
-         enddo
-      endif
 
 c FH modif sur le cdrag temperature
@@ -737,4 +761,5 @@
 #else
        DO j=1, knon
+         i = ni(j) 
          t2m(i,nsrf)=0.
          q2m(i,nsrf)=0.
@@ -1219,4 +1244,7 @@
       END
       SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
+cIM 261103
+     .                  ksta, ksta_ter,
+cIM 261103
      .                  ts, rugos,
      .                  u,v,t,q,
@@ -1279,5 +1307,13 @@
       PARAMETER (prandtl=0.4)
       REAL kstable ! diffusion minimale (situation stable)
-      PARAMETER (kstable=1.0e-10)
+      ! GKtest
+      ! PARAMETER (kstable=1.0e-10)
+      REAL ksta, ksta_ter
+cIM: 261103     REAL kstable_ter, kstable_sinon
+cIM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)
+cIM: 261103     PARAMETER (kstable_ter = 1.0e-8)
+cIM: 261103   PARAMETER (kstable_ter = 1.0e-10)
+cIM: 261103   PARAMETER (kstable_sinon = 1.0e-10)
+      ! fin GKtest
       REAL mixlen ! constante controlant longueur de melange
       PARAMETER (mixlen=35.0)
@@ -1371,4 +1407,13 @@
          gamt(2) = -2.5E-03
       ENDIF
+cIM cf JLD/ GKtest
+      IF ( nsrf .EQ. is_ter ) THEN
+cIM 261103     kstable = kstable_ter
+        kstable = ksta_ter
+      ELSE
+cIM 261103     kstable = kstable_sinon
+        kstable = ksta
+      ENDIF
+cIM cf JLD/ GKtest fin
 c
 c Calculer les geopotentiels de chaque couche
@@ -1614,7 +1659,12 @@
       DO k = 2, klev
       DO i = 1, knon
-      IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
-     .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
-     .     (zdthmin(i).GT.seuil) )THEN ! si l'inversion est trop faible
+cIM cf FH/GK   IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
+cIM cf FH/GK  .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
+      !IM cf JLD/ GKtest TERkz2
+      ! IF ( (nsrf.EQ.is_ter) .OR.  ! si on est sur la terre
+      ! fin GKtest
+      IF ( (nsrf.EQ.is_oce) .AND.  ! si on est sur ocean et si 
+     .     ( (invb(i).EQ.klev) .OR.      ! s'il n'y a pas d'inversion
+     .     (zdthmin(i).GT.seuil) ) )THEN ! si l'inversion est trop faible
          zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
      .                       /(paprs(i,2)-paprs(i,klev+1)) ))**2
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 486)
@@ -442,4 +442,33 @@
   call getin('cdhmax',cdhmax)
 
+!261103
+!
+!Config Key  = ksta
+!Config Desc =
+!Config Def  = 1.0e-10
+!Config Help =
+!
+  ksta = 1.0e-10
+  call getin('ksta',ksta)
+
+!
+!Config Key  = ksta_ter
+!Config Desc =
+!Config Def  = 1.0e-10
+!Config Help =
+!
+  ksta_ter = 1.0e-10
+  call getin('ksta_ter',ksta_ter)
+
+!
+!Config Key  = ok_kzmin
+!Config Desc =
+!Config Def  = .true.
+!Config Help =
+!
+  ok_kzmin = .true.
+  call getin('ok_kzmin',ok_kzmin)
+
+!261103
 !
 !Config Key  = 
@@ -495,4 +524,9 @@
   write(numout,*)' top_height = ',top_height 
   write(numout,*)' overlap = ',overlap 
+  write(numout,*)' cdmmax = ',cdmmax 
+  write(numout,*)' cdhmax = ',cdhmax 
+  write(numout,*)' ksta = ',ksta 
+  write(numout,*)' ksta_ter = ',ksta_ter 
+  write(numout,*)' ok_kzmin = ',ok_kzmin 
 
   return
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F	(revision 486)
@@ -1,2 +1,5 @@
+c
+c $Header$
+c
       SUBROUTINE cv3_param(nd,delt)
       implicit none
@@ -1459,4 +1462,5 @@
       real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)
       real wgh
+      real zm(nloc,na)
       logical lwork(nloc)
 
@@ -1498,4 +1502,5 @@
        enddo
       enddo
+      zm(:,:)=0.
 
 c=====================================================================
@@ -1579,5 +1584,6 @@
       vent(il,i,i)=v(il,nk(il))
       elij(il,i,i)=clw(il,i)
-      sij(il,i,i)=1.0
+cMAF      sij(il,i,i)=1.0
+      sij(il,i,i)=0.0
       end if
  740  continue
@@ -1753,5 +1759,6 @@
         vent(il,i,i)=v(il,nk(il))
         elij(il,i,i)=clw(il,i)
-        sij(il,i,i)=1.0
+cMAF        sij(il,i,i)=1.0
+        sij(il,i,i)=0.0
        endif
       enddo ! il
@@ -1767,5 +1774,24 @@
 
 789   continue
-      
+c      
+c MAF: renormalisation de MENT
+      do jm=1,nd
+        do im=1,nd
+          do il=1,ncum
+          zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
+         end do
+        end do
+      end do
+c
+      do jm=1,nd
+        do im=1,nd
+          do il=1,ncum
+          if(zm(il,im).ne.0.) then
+          ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
+          endif
+         end do
+       end do
+      end do
+c
       do jm=1,nd
        do im=1,nd
@@ -2079,5 +2105,6 @@
       do j=1,ntra
       trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
-     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
+ctestmaf     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
+     :            +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
       trap(il,i,j)=trap(il,i,j)/mp(il,i)
       end do
@@ -2636,9 +2663,9 @@
           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
      :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
-     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
+     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
          else
           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
      :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
-     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
+     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
          endif
         endif ! i
@@ -2692,5 +2719,5 @@
         ex=0.1*ment(il,inb(il),inb(il)) 
      :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
-     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
+     :      /(ph(il,inb(il))-ph(il,inb(il)+1))
         ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
         ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
@@ -2792,8 +2819,9 @@
       enddo
 
-      do i=1,nl
+      do i=2,nl
        do k=i,nl
         do il=1,ncum
-         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
+ctest         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
+         if (i.le.inb(il).and.k.le.inb(il)) then
             upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
             dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
@@ -3024,5 +3052,4 @@
 
         do 2100 j=1,ntra
-c oct3         do 2110 k=1,nl
          do 2110 k=1,nd ! oct3
           do 2120 i=1,ncum
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/hydrol.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/hydrol.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/hydrol.F	(revision 486)
@@ -1,2 +1,5 @@
+c
+c $Header$
+c
       SUBROUTINE hydrol(dtime,pctsrf,rain_fall,snow_fall,evap,
      .                  agesno, tsol,qsol,snow,runoff)
@@ -76,5 +79,7 @@
 c je limite la temperature a RTT-1.8 (il faudrait aussi prendre l'eau de
 c la fonte) (Laurent Li, le 14mars98):
-      tsol(i,is) = MIN(tsol(i,is),RTT-1.8)
+cIM cf GK   tsol(i,is) = MIN(tsol(i,is),RTT-1.8)
+cIM cf GK : la glace fond a 0C, non pas a -1.8
+            tsol(i,is) = MIN(tsol(i,is),RTT)
 c
 ccc         ELSE
@@ -95,5 +100,7 @@
 c je limite la temperature a RTT-1.8 (il faudrait aussi prendre l'eau de
 c la fonte) (Laurent Li, le 14mars98):
-      tsol(i,is) = MIN(tsol(i,is),RTT-1.8)
+cIM cf GK   tsol(i,is) = MIN(tsol(i,is),RTT-1.8)
+cIM cf GK : la glace fond a 0C, non pas a -1.8
+            tsol(i,is) = MIN(tsol(i,is),RTT)
 c
 ccc         ELSE
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histday.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histday.h	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histday.h	(revision 486)
@@ -233,23 +233,4 @@
      .                capemaxcels, zsto,zout)
 
-         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)
 
 c=================================================================
@@ -438,22 +419,4 @@
      .                "ave(X)", zsto,zout)
 c
-cccIM   
-         CALL histdef(nid_day, "SWupTOA", "SWup at TOA","W/m2",
-     .                iim,jjmp1,nhori, 1,1,1,-99,
-     .                32, "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_day, "SWupSFC", "SWup at surface","W/m2",
-     .                iim,jjmp1,nhori, 1,1,1,-99,
-     .                32, "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_day, "SWdnTOA", "SWdn at TOA","W/m2",
-     .                iim,jjmp1,nhori, 1,1,1,-99,
-     .                32, "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_day, "SWdnSFC", "SWdn at surface","W/m2",
-     .                iim,jjmp1,nhori, 1,1,1,-99,
-     .                32, "ave(X)", zsto,zout)
-c
-         CALL histend(nid_day)
 c
 c=================================================================
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/ini_histmth.h	(revision 486)
@@ -116,6 +116,10 @@
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
      .                "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_mth, "solldown", "Down. IR rad. at surface", 
+cIM: 071003
+         CALL histdef(nid_mth,"LWdnSFC","Down. IR rad. at surface", 
+     .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+cIM: 071003
+         CALL histdef(nid_mth,"LWupSFC","Upwd. IR rad. at surface", 
      .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32, 
      .                "ave(X)", zsto,zout)
@@ -263,6 +267,22 @@
      .                "ave(X)", zsto,zout)
 c
-         CALL histdef(nid_mth, "cldq", "Cloud liquid water path", "-",
-     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+         CALL histdef(nid_mth,"cldq","Cloud liquid water path","Kg/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+cIM: 071003
+         CALL histdef(nid_mth,"lwp","Cloud water path","Kg/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth,"iwp","Cloud ice water path","Kg/m2",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+cIM: 071003
+         CALL histdef(nid_mth,"lwcon","Cloud water content","Kg/Kg",
+     .                iim,jjmp1,nhori, klev,1,klev, nvert, 32, 
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth,"iwcon","Cloud ice water content","Kg/Kg",
+     .                iim,jjmp1,nhori, klev,1,klev, nvert, 32, 
      .                "ave(X)", zsto,zout)
 c
@@ -334,4 +354,8 @@
      .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
      .                "ave(X)", zsto,zout)
+cIM: 071003
+         CALL histdef(nid_mth,"wvap","Water vapor mixing ratio","Kg/Kg",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
 c
          CALL histdef(nid_mth, "geop", "Geopotential height", "m",
@@ -374,4 +398,8 @@
 c
          CALL histdef(nid_mth, "oliq", "Liquid water content", "kg/kg",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_mth, "dtphy", "Physics dT", "K/s",
      .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
      .                "ave(X)", zsto,zout)
@@ -400,4 +428,8 @@
      .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
      .                "ave(X)", zsto,zout)
+cIM: 071003
+         CALL histdef(nid_mth, "dtlschr",
+     $       "Large-scale condensational heating rate", "K/s",iim,jjmp1
+     $       ,nhori, klev,1,klev,nvert, 32,"ave(X)", zsto,zout)
 c
          CALL histdef(nid_mth, "dqlsc", "Condensation dQ", "Kg/Kg/s",
@@ -442,5 +474,5 @@
      .                "ave(X)", zsto,zout)
 c
-         CALL histdef(nid_mth, "dtsw0", "SW radiation dT", "K/s",
+         CALL histdef(nid_mth, "dtsw0", "CS SW radiation dT", "K/s",
      .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
      .                "ave(X)", zsto,zout)
@@ -450,5 +482,5 @@
      .                "ave(X)", zsto,zout)
 c
-         CALL histdef(nid_mth, "dtlw0", "LW radiation dT", "K/s",
+         CALL histdef(nid_mth,"dtlw0","CS LW radiation dT","K/s",
      .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
      .                "ave(X)", zsto,zout)
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 486)
@@ -694,5 +694,12 @@
      alb_new(1 : knon)  = alb_neig(1 : knon)*zfra(1:knon) + &
     &                     0.6 * (1.0-zfra(1:knon))
-!!     alb_new(1 : knon)  = 0.6
+!IM cf FH/GK     alb_new(1 : knon)  = 0.6
+!       alb_new(1 : knon)  = 0.82
+!IM cf JLD/ GK
+!IM: 211003 Ksta0.77      alb_new(1 : knon)  = 0.77
+!IM: KstaTER0.8 & LMD_ARMIP5    alb_new(1 : knon)  = 0.8
+!IM: KstaTER0.77 & LMD_ARMIP6    
+        alb_new(1 : knon)  = 0.77
+
 !
 ! Rugosite
@@ -1144,5 +1151,6 @@
 
 !IM cf. JP +++
-    albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2.
+!IM BUG BUG BUG 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 ---
 
@@ -2599,5 +2607,7 @@
 ! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige
 ! en exces "s'ecoule" (calving)
-  real, parameter :: snow_max=1.
+!  real, parameter :: snow_max=1.
+!IM cf JLD/GK
+  real, parameter :: snow_max=3000.
   integer :: i
   real, dimension(klon) :: zx_mh, zx_nh, zx_oh
@@ -2614,4 +2624,7 @@
 ! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15)
   REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
+!IM cf JLD/ GKtest
+  REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)
+! fin GKtest
 !
   logical, save         :: check = .FALSE.
@@ -2702,5 +2715,13 @@
       tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno  
 !IM cf JLD OK     
-      IF (nisurf == is_sic .OR. nisurf == is_lic ) tsurf_new(i) = RTT
+!IM cf JLD/ GKtest fonte aussi pour la glace
+!      IF (nisurf == is_sic .OR. nisurf == is_lic ) tsurf_new(i) = RTT
+      IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
+        fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)
+        ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
+        bil_eau_s(i) = bil_eau_s(i) + fq_fonte
+        tsurf_new(i) = RTT
+      ENDIF
+! fin GKtest
       d_ts(i) = tsurf_new(i) - tsurf(i)
 !      zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
@@ -2729,6 +2750,7 @@
       run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
       qsol(i) = MIN(qsol(i), max_eau_sol) 
-    else
-      run_off(i) = run_off(i) + MAX(bil_eau_s(i), 0.0)
+!IM : 0601003 else
+!IM: run_off(i)  
+!IM : 061003   run_off(i) = run_off(i) + MAX(bil_eau_s(i), 0.0)
     endif
   enddo
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/isccp_cloud_types.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/isccp_cloud_types.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/isccp_cloud_types.F	(revision 486)
@@ -281,8 +281,8 @@
 c          write(6,'(a10)') 'ncol='
 c          write(6,'(8I10)') ncol
-           write(6,'(a10)') 'top_height='
-           write(6,'(8I10)') top_height
-           write(6,'(a10)') 'overlap='
-           write(6,'(8I10)') overlap
+c          write(6,'(a10)') 'top_height='
+c          write(6,'(8I10)') top_height
+c          write(6,'(a10)') 'overlap='
+c          write(6,'(8I10)') overlap
 c          write(6,'(a10)') 'emsfc_lw='
 c          write(6,'(8f10.2)') emsfc_lw
@@ -656,7 +656,5 @@
     
           ! Reset threshold 
-
-cIM pas besoin..memes val. de ran!      
-c         call ran0_vec(npoints,seed,ran)
+          call ran0_vec(npoints,seed,ran)
 	   
           do j=1,npoints
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/newmicro.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/newmicro.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/newmicro.F	(revision 486)
@@ -1,5 +1,8 @@
       SUBROUTINE newmicro (paprs, pplay,ok_newmicro,
      .                  t, pqlwp, pclc, pcltau, pclemi,
-     .                  pch, pcl, pcm, pct, pctlwp)
+cIM    .                  pch, pcl, pcm, pct, pctlwp)
+     .                  pch, pcl, pcm, pct, pctlwp,
+     .                  xflwp, xfiwp, xflwc, xfiwc)
+
       IMPLICIT none
 c======================================================================
@@ -36,5 +39,9 @@
 C
       INTEGER i, k
-      REAL zflwp, zradef, zfice, zmsac
+cIM: 091003   REAL zflwp, zradef, zfice, zmsac
+      REAL zflwp(klon), zradef, zfice, zmsac
+cIM: 091003 rajout
+      REAL xflwp(klon), xfiwp(klon)
+      REAL xflwc(klon,klev), xfiwc(klon,klev)
 c
       REAL radius, rad_chaud
@@ -53,5 +60,6 @@
       logical ok_newmicro
 c     parameter (ok_newmicro=.FALSE.)
-      real rel, tc, rei, zfiwp
+cIM: 091003   real rel, tc, rei, zfiwp
+      real rel, tc, rei, zfiwp(klon)
       real k_liq, k_ice0, k_ice, DF
       parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g
@@ -62,10 +70,17 @@
 c Calculer l'epaisseur optique et l'emmissivite des nuages
 c
+cIM inversion des DO
+      DO i = 1, klon
+       xflwp(i)=0.
+       xfiwp(i)=0.
       DO k = 1, klev
-      DO i = 1, klon
+c
+       xflwc(i,k)=0.
+       xfiwc(i,k)=0.
+c
          rad_chaud = rad_chau1
          IF (k.LE.3) rad_chaud = rad_chau2
          pclc(i,k) = MAX(pclc(i,k), seuil_neb)
-         zflwp = 1000.*pqlwp(i,k)/RG/pclc(i,k)
+         zflwp(i) = 1000.*pqlwp(i,k)/RG/pclc(i,k)
      .          *(paprs(i,k)-paprs(i,k+1))
          zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
@@ -74,6 +89,6 @@
          radius = rad_chaud * (1.-zfice) + rad_froid * zfice
          coef = coef_chau * (1.-zfice) + coef_froi * zfice
-         pcltau(i,k) = 3.0/2.0 * zflwp / radius
-         pclemi(i,k) = 1.0 - EXP( - coef * zflwp)
+         pcltau(i,k) = 3.0/2.0 * zflwp(i) / radius
+         pclemi(i,k) = 1.0 - EXP( - coef * zflwp(i))
 
          if (ok_newmicro) then
@@ -84,8 +99,20 @@
          zfice = MIN(MAX(zfice,0.0),1.0)
 
-         zflwp = 1000.*(1.-zfice)*pqlwp(i,k)/pclc(i,k)
-     :          *(paprs(i,k)-paprs(i,k+1))/RG
-         zfiwp = 1000.*zfice*pqlwp(i,k)/pclc(i,k)
-     :          *(paprs(i,k)-paprs(i,k+1))/RG
+         zflwp(i) = 1000.*(1.-zfice)*pqlwp(i,k)/pclc(i,k)
+     :          *(paprs(i,k)-paprs(i,k+1))/RG
+         zfiwp(i) = 1000.*zfice*pqlwp(i,k)/pclc(i,k)
+     :          *(paprs(i,k)-paprs(i,k+1))/RG
+
+         xflwp(i) = xflwp(i)+ (1.-zfice)*pqlwp(i,k)
+     :          *(paprs(i,k)-paprs(i,k+1))/RG
+         xfiwp(i) = xfiwp(i)+ zfice*pqlwp(i,k)
+     :          *(paprs(i,k)-paprs(i,k+1))/RG
+
+cIM Total Liquid/Ice water content
+         xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)
+         xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)
+cIM In-Cloud Liquid/Ice water content
+c        xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
+c        xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
 
 c -- effective cloud droplet radius (microns):
@@ -107,8 +134,8 @@
 c  for ice clouds, Ebert & Curry (1992)] 
 
-         if (zflwp.eq.0.) rel = 1. 
-         if (zfiwp.eq.0. .or. rei.le.0.) rei = 1. 
-         pcltau(i,k) = 3.0/2.0 * ( zflwp/rel )
-     .             + zfiwp * (3.448e-03  + 2.431/rei)
+         if (zflwp(i).eq.0.) rel = 1. 
+         if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1. 
+         pcltau(i,k) = 3.0/2.0 * ( zflwp(i)/rel )
+     .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
 
 c -- cloud infrared emissivity:
@@ -121,5 +148,5 @@
 
          pclemi(i,k) = 1.0
-     .      - EXP( - coef_chau*zflwp - DF*k_ice*zfiwp )
+     .      - EXP( - coef_chau*zflwp(i) - DF*k_ice*zfiwp(i) )
 
          endif ! ok_newmicro
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 486)
@@ -40,5 +40,6 @@
 #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
@@ -554,4 +555,19 @@
       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'
@@ -574,4 +590,54 @@
       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/physiq.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 486)
@@ -132,4 +132,8 @@
 c      PARAMETER (ok_mensuel=.true.)
 c
+      LOGICAL ok_mensuelNMC ! sortir le fichier mensuel niveaux NMC
+      PARAMETER (ok_mensuelNMC=.true.)
+c     save ok_mensuelNMC
+c
       LOGICAL ok_instan ! sortir le fichier instantane
       save ok_instan
@@ -187,23 +191,33 @@
       REAL d_ps(klon)
 
-cccIM
-      INTEGER klevp1
-      PARAMETER(klevp1=klev+1)
+      INTEGER klevp1, klevm1
+      PARAMETER(klevp1=klev+1,klevm1=klev-1)
 #include "raddim.h"
-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
+c vents meridien et zonal a un niveau de pression
+      real u1000(klon), v1000(klon) !vents a 1000 hPa
+      real u925(klon), v925(klon)   !vents a  925 hPa
+      real u850(klon),v850(klon)    !vents a  850 hPa
+      real u700(klon),v700(klon)
+      real u600(klon),v600(klon)
+      real u500(klon),v500(klon)
+      real u400(klon),v400(klon)
+      real u300(klon),v300(klon)
+      real u250(klon),v250(klon)
+      real u200(klon),v200(klon)
+      real u150(klon),v150(klon)
+      real u100(klon),v100(klon)
+      real u70(klon),v70(klon)
+      real u50(klon),v50(klon)
+      real u30(klon),v30(klon)
+      real u20(klon),v20(klon)
+      real u10(klon),v10(klon)
+      real phi500(klon),w500(klon)
+c prw: precipitable water
       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
@@ -214,20 +228,40 @@
       REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
 
-      INTEGER kinv, linv
-
-cIM ISCCP simulator BEGIN
-      INTEGER igfi2D(iim,jjmp1)
+      INTEGER linv, kp1
+c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
+c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
+      REAL flwp(klon), fiwp(klon)
+      REAL flwc(klon,klev), fiwc(klon,klev)
+      REAL flwp_c(klon), fiwp_c(klon)
+      REAL flwc_c(klon,klev), fiwc_c(klon,klev)
+      REAL flwp_s(klon), fiwp_s(klon)
+      REAL flwc_s(klon,klev), fiwc_s(klon,klev)
+
+c ISCCP simulator v3.4
+c dans clesphys.h top_height, overlap
 cv3.4
       INTEGER debug, debugcol
       INTEGER npoints
       PARAMETER(npoints=klon) 
-      INTEGER sunlit(klon)
-
+c
+      INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
+      INTEGER nregISCtot
+      PARAMETER(nregISCtot=1) 
+c
+c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
+c y compris pour 1 point
+c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
+c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
+      INTEGER imin_debut, nbpti
+      INTEGER jmin_debut, nbptj 
+c
+      REAL nbsunlit(nregISCtot,klon)  !nbsunlit : moyenne de sunlit
       INTEGER ncol, seed(klon)
 
-cIM dans clesphys.h top_height, overlap
+c ncol = nb. de sous-colonnes pour chaque maille du GCM 
 c     PARAMETER(ncol=100)
 c     PARAMETER(ncol=625)
-      PARAMETER(ncol=10)
+c     PARAMETER(ncol=10)
+      PARAMETER(ncol=25)
       REAL tautab(0:255)
       INTEGER invtau(-20:45000)
@@ -235,13 +269,15 @@
       PARAMETER(emsfc_lw=0.99)
       REAL    ran0                      ! type for random number fuction
-
+c
+      REAL cldtot(klon,klev)
+c variables de haut en bas pour le simulateur ISCCP
+      REAL dtau_s(klon,klev) !tau nuages startiformes
+      REAL dtau_c(klon,klev) !tau nuages convectifs
+      REAL dem_s(klon,klev)  !emissivite nuages startiformes 
+      REAL dem_c(klon,klev)  !emissivite nuages convectifs
+c
+c variables de haut en bas pour le simulateur ISCCP
       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)
@@ -253,5 +289,5 @@
       REAL dem_cH2B(klon,klev)
 
-c output from ISCCP
+c output from ISCCP simulator
       REAL fq_isccp(klon,7,7)
       REAL totalcldarea(klon) 
@@ -260,23 +296,14 @@
       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)
+c
+      INTEGER l, ni, nj, kmax, lmax
       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
+      INTEGER iimx7, jjmx7, jjmp1x7
+      PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1, 
+     .jjmp1x7=jjmp1*lmaxm1)
+      REAL fq4d(iim,jjmp1,kmaxm1,lmaxm1)
+      REAL fq3d(iimx7, jjmp1x7)
 c
       INTEGER iw, iwmax
@@ -285,22 +312,16 @@
       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) 
+c
+cIM: nbregdyn = nbre regions pour calculs statistiques sur output du ISCCP
+cIM: dynamiques  
+      INTEGER nreg, nbregdyn
+      PARAMETER(nbregdyn=5)
+      REAL histoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
+      REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
+      REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbregdyn) 
       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)
+      INTEGER pct_ocean(klon,nbregdyn)
       REAL rlonPOS(klon) 
-c     CHARACTER*4 pdirect
  
 c sorties ISCCP
@@ -321,29 +342,40 @@
 #endif
 
+c sorties statistiques regime dynamique
+      logical ok_regdyn
+      real ecrit_regdyn
+      integer nid_regdyn
+      save ok_regdyn, ecrit_regdyn, nid_regdyn
+
+#undef histREGDYN
+#define histREGDYN
+#ifdef histREGDYN
+c     data ok_regdyn,ecrit_regdyn/.true.,0.125/
+c     data ok_regdyn,ecrit_regdyn/.true.,1./
+       data ok_regdyn/.true./
+#else
+      data ok_regdyn/.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./
+      DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
+      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
+
+c cldtopres pression au sommet des nuages
+      REAL cldtopres(lmaxm1)
+      DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
 
       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
-
+c taulev: numero du niveau de tau dans les sorties ISCCP
+      CHARACTER *4 taulev(kmaxm1)
+      DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
+
+      REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
+      INTEGER nhorix7
+cIM: region='3d' <==> sorties en global
+      CHARACTER*3 region
+      PARAMETER(region='3d')
+c
       logical ok_hf
       real ecrit_hf
@@ -513,7 +545,6 @@
       REAL yu1(klon)            ! vents dans la premiere couche U
       REAL yv1(klon)            ! vents dans la premiere couche V
-cIM cf JLD
-      REAL ffonte(klon,nbsrf)         !Flux thermique utilise pour fondre la neige
-      REAL fqcalving(klon,nbsrf)      !Flux d'eau "perdue" par la surface 
+      REAL ffonte(klon,nbsrf)    !Flux thermique utilise pour fondre la neige
+      REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface 
 c                               !et necessaire pour limiter la
 c                               !hauteur de neige, en kg/m2/s
@@ -539,5 +570,4 @@
       REAL dlw(klon)    ! derivee infra rouge
       REAL bils(klon) ! bilan de chaleur au sol
-cIM cf. JLD
       REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
 C                   type de sous-surface et pondere par la fraction
@@ -574,5 +604,4 @@
       EXTERNAL angle     ! calculer angle zenithal du soleil
       EXTERNAL alboc     ! calculer l'albedo sur ocean
-      EXTERNAL albsno    ! calculer albedo sur neige
       EXTERNAL ajsec     ! ajustement sec
       EXTERNAL clmain    ! couche limite 
@@ -601,4 +630,6 @@
       EXTERNAL ecrirega  ! ecrire le fichier binaire regional
       EXTERNAL ecriregs  ! ecrire le fichier binaire regional
+cIM
+      EXTERNAL haut2bas  !variables de haut en bas
 c
 c Variables locales
@@ -685,6 +716,5 @@
       REAL cape(klon)           ! CAPE
       SAVE cape
-cccIM
-      CHARACTER*40 capemaxcels
+      CHARACTER*40 capemaxcels  !max(CAPE)
 
       REAL pbase(klon)          ! cloud base pressure
@@ -739,4 +769,5 @@
       REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
       REAL d_t_lif(klon,klev)
+      REAL d_u_oli(klon,klev), d_v_oli(klon,klev) !tendances dues a oro et lif 
 
       REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
@@ -792,10 +823,11 @@
 c
       INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
-      REAL zx_tmp_fi2d(klon)
+      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
+      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
       REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
       REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
 c
-      INTEGER nid_day, nid_mth, nid_ins
-      SAVE nid_day, nid_mth, nid_ins
+      INTEGER nid_day, nid_mth, nid_ins, nid_nmc
+      SAVE nid_day, nid_mth, nid_ins, nid_nmc
 c
       INTEGER nhori, nvert
@@ -841,10 +873,10 @@
       REAL ZRCPD
 c-jld ec_conser
-cIM
-      REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
-      REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
-      REAL zt2m(klon), zq2m(klon)
-      REAL zu10m(klon), zv10m(klon)
-      CHARACTER*40 t2mincels, t2maxcels 
+cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
+      REAL t2m(klon,nbsrf), q2m(klon,nbsrf)   !temperature, humidite a 2m
+      REAL u10m(klon,nbsrf), v10m(klon,nbsrf) !vents a 10m
+      REAL zt2m(klon), zq2m(klon)             !temp., hum. 2m moyenne s/ 1 maille
+      REAL zu10m(klon), zv10m(klon)           !vents a 10m moyennes s/1 maille
+      CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
 c
 c Declaration des constantes et des fonctions thermodynamiques
@@ -1031,16 +1063,26 @@
 c   Initialisation des sorties
 c=============================================================
+#ifdef histhf
+#include "ini_histhf.h"
+#endif
+
+#include "ini_histday.h"
+#include "ini_histmth.h"
+
+#undef histmthNMC
+#define histmthNMC
+#ifdef histmthNMC
+#include "ini_histmthNMC.h"
+#endif
+
+#include "ini_histins.h"
+
+#ifdef histREGDYN
+#include "ini_histREGDYN.h"
+#endif
 
 #ifdef histISCCP
 #include "ini_histISCCP.h"
 #endif
-
-#ifdef histhf
-#include "ini_histhf.h"
-#endif
-
-#include "ini_histday.h"
-#include "ini_histmth.h"
-#include "ini_histins.h"
 
 cXXXPB Positionner date0 pour initialisation de ORCHIDEE
@@ -1253,8 +1295,5 @@
        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
+       nbsunlit(1,i)=FLOAT(sunlit(i))
       ENDDO
 cIM END
@@ -1289,5 +1328,6 @@
      e            julien, rmu0, co2_ppm, 
      e            ok_veget, ocean, npas, nexca, ftsol,
-     $            soil_model,cdmmax, cdhmax, ftsoil, qsol, 
+     $            soil_model,cdmmax, cdhmax,
+     $            ksta, ksta_ter, ok_kzmin, ftsoil, qsol, 
      $            paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
      $            fluxlat,
@@ -1789,55 +1829,33 @@
       enddo
 
-cIM ISCCP simulator BEGIN
+cIM calcul nuages par le simulateur ISCCP
       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)
-
+     .            cldh_c, cldl_c, cldm_c, cldt_c, cldq_c,
+     .            flwp_c, fiwp_c, flwc_c, fiwc_c)
+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     
+     .            cldh_s, cldl_s, cldm_s, cldt_s, cldq_s,
+     .            flwp_s, fiwp_s, flwc_s, fiwc_s)
+c
       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 
+      CALL haut2bas(klon, klev, pplay, pfull)
+      CALL haut2bas(klon, klev, q_seri, qv)
+      CALL haut2bas(klon, klev, cldtot, cc)
+      CALL haut2bas(klon, klev, rnebcon, conv)
+      CALL haut2bas(klon, klev, dtau_s, dtau_sH2B)
+      CALL haut2bas(klon, klev, dtau_c, dtau_cH2B)
+      CALL haut2bas(klon, klev, t_seri, at)
+      CALL haut2bas(klon, klev, dem_s, dem_sH2B)
+      CALL haut2bas(klon, klev, dem_c, dem_cH2B)
+      CALL haut2bas(klon, klevp1, paprs, phalf)
 
 c     open(99,file='tautab.bin',access='sequential',
@@ -1855,6 +1873,9 @@
       close(99)
 c
+cIM: calcul coordonnees regions pour statistiques distribution 
+cIM: nuages en ftion du regime dynamique pour regions oceaniques
+       IF (ok_regdyn) THEN !histREGDYN
        nsrf=3
-       DO nreg=1, nbreg
+       DO nreg=1, nbregdyn
        DO i=1, klon
 
@@ -1867,8 +1888,5 @@
 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
@@ -1879,5 +1897,4 @@
 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
@@ -1887,5 +1904,4 @@
            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
@@ -1895,5 +1911,4 @@
           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
@@ -1903,5 +1918,4 @@
          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
@@ -1911,9 +1925,8 @@
          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
+        ENDIF !nbregdyn
 c TROP
 c        IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
@@ -1924,143 +1937,43 @@
 
         ENDIF !pctsrf
-c      ENDDO 
        ENDDO !klon
-       ENDDO !nbreg
+       ENDDO !nbregdyn
+       ENDIF !ok_regdyn
  
 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
+      DO nreg = 1, nbregdyn
+       DO k = 1, kmaxm1
+        DO l = 1, lmaxm1
+         DO iw = 1, iwmax
+          nhistoWt(k,l,iw,nreg)=0.
+         ENDDO !iw
+        ENDDO !l
+       ENDDO !k
+      ENDDO !nreg
 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
+c
+cIM: initialisation de seed
         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
+      ENDIF !debut
+cIM: pas de debug, debugcol
       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) 
+c o500 ==> distribution nuage ftion du regime dynamique a 500 hPa
+        DO k=1, klevm1
+        kp1=k+1
+c       PRINT*,'k, presnivs',k,presnivs(k), presnivs(kp1)
+        if(presnivs(k).GT.50000.AND.presnivs(kp1).LT.50000.) THEN
+         DO i=1, klon
+          o500(i)=omega(i,k)*RDAY/100.
+c         if(i.EQ.1) print*,' 500hPa lev',k,presnivs(k),presnivs(kp1)
+         ENDDO
+         GOTO 1000
+        endif 
+1000  continue
+      ENDDO
+
       CALL ISCCP_CLOUD_TYPES(
      &     debug,
@@ -2073,13 +1986,5 @@
      &     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,
@@ -2088,13 +1993,5 @@
      &     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,
@@ -2103,259 +2000,56 @@
      &     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 l=1, lmaxm1
+       DO k=1, kmaxm1
         DO i=1, iim
-          fq4d(i,1,k,l)=fq_isccp(1,k,l) 
-cc         PRINT*,'first j=1 i =',i
+         fq4d(i,1,k,l)=fq_isccp(1,k,l) 
         ENDDO
         DO j=2, jjm
-          DO i=1, iim
-cERROR ??         ig=i+iim*(j-1)
+         DO i=1, iim
           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
+         fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l) 
         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  
+c
+      DO l=1, lmaxm1
+       DO k=1, kmaxm1  
         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)
-         
+           ni=(i-1)*lmaxm1+l
+           nj=(j-1)*kmaxm1+k
+           fq3d(ni,nj)=fq4d(i,j,k,l)
          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
+      ENDDO
+
 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
 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,
+      CALL histo_o500_pctau(nbregdyn,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
-
-
+c nhistoWt = somme de toutes les nhistoW
+      DO nreg=1, nbregdyn
+       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)
+         ENDDO
+        ENDDO
+       ENDDO
+      ENDDO
+c
       ENDIF !ok_isccp
-cIM ISCCP simulator END
 
 c   On prend la somme des fractions nuageuses et des contenus en eau
@@ -2363,6 +2057,6 @@
       cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
 
-
       ENDIF
+
 c
 c 2. NUAGES STARTIFORMES
@@ -2423,5 +2117,6 @@
       CALL newmicro (paprs, pplay,ok_newmicro,
      .            t_seri, cldliq, cldfra, cldtau, cldemi,
-     .            cldh, cldl, cldm, cldt, cldq)
+     .            cldh, cldl, cldm, cldt, cldq,
+     .            flwp, fiwp, flwc, fiwc)
       else
       CALL nuage (paprs, pplay,
@@ -2450,5 +2145,4 @@
 !      albsollw = albsollw1
       CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
-cIM  e            (dist, rmu0, fract, co2_ppm, solaire,
      e            (dist, rmu0, fract, 
      e             paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
@@ -2458,5 +2152,4 @@
      s             topsw,toplw,solsw,sollw,
      s             sollwdown,
-cccIMs             topsw0,toplw0,solsw0,sollw0)
      s             topsw0,toplw0,solsw0,sollw0,
      s             swdn0, swdn, swup0, swup     )
@@ -2658,4 +2351,5 @@
      s                   ve, vq, ue, uq)
 c
+c
 c Accumuler les variables a stocker dans les fichiers histoire:
 c
@@ -2692,5 +2386,4 @@
       END IF 
 C
-cccIM cf. FH
 c=======================================================================
 c   SORTIES
@@ -2699,33 +2392,85 @@
 c   Interpollation sur quelques niveaux de pression
 c   -----------------------------------------------
-
+c
+cIM sorties sur les 17 niveaux de pression du NMC
+c 1000 hPa
+      call plevel(klon,klev,.true. ,pplay,100000.,u_seri,u1000)
+      call plevel(klon,klev,.false.,pplay,100000.,v_seri,v1000)
+c 925 hPa
+      call plevel(klon,klev,.true. ,pplay,92500.,u_seri,u925)
+      call plevel(klon,klev,.false.,pplay,92500.,v_seri,v925)
+c 850 hPa
       call plevel(klon,klev,.true. ,pplay,85000.,u_seri,u850)
       call plevel(klon,klev,.false.,pplay,85000.,v_seri,v850)
+c 700 hPa
+      call plevel(klon,klev,.true. ,pplay,70000.,u_seri,u700)
+      call plevel(klon,klev,.false.,pplay,70000.,v_seri,v700)
+c 600 hPa
+      call plevel(klon,klev,.true. ,pplay,60000.,u_seri,u600)
+      call plevel(klon,klev,.false.,pplay,60000.,v_seri,v600)
+c 500 hPa
       call plevel(klon,klev,.true. ,pplay,50000.,u_seri,u500)
       call plevel(klon,klev,.false.,pplay,50000.,v_seri,v500)
+c 400 hPa
+      call plevel(klon,klev,.true. ,pplay,40000.,u_seri,u400)
+      call plevel(klon,klev,.false.,pplay,40000.,v_seri,v400)
+c 300 hPa
+      call plevel(klon,klev,.true. ,pplay,30000.,u_seri,u300)
+      call plevel(klon,klev,.false.,pplay,30000.,v_seri,v300)
+c 250 hPa
+      call plevel(klon,klev,.true. ,pplay,25000.,u_seri,u250)
+      call plevel(klon,klev,.false.,pplay,25000.,v_seri,v250)
+c 200 hPa
       call plevel(klon,klev,.true. ,pplay,20000.,u_seri,u200)
       call plevel(klon,klev,.false.,pplay,20000.,v_seri,v200)
+c 150 hPa
+      call plevel(klon,klev,.true. ,pplay,15000.,u_seri,u150)
+      call plevel(klon,klev,.false.,pplay,15000.,v_seri,v150)
+c 100 hPa
+      call plevel(klon,klev,.true. ,pplay,10000.,u_seri,u100)
+      call plevel(klon,klev,.false.,pplay,10000.,v_seri,v100)
+c 70 hPa
+      call plevel(klon,klev,.true. ,pplay,7000.,u_seri,u70)
+      call plevel(klon,klev,.false.,pplay,7000.,v_seri,v70)
+c 50 hPa
+      call plevel(klon,klev,.true. ,pplay,5000.,u_seri,u50)
+      call plevel(klon,klev,.false.,pplay,5000.,v_seri,v50)
+c 30 hPa
+      call plevel(klon,klev,.true. ,pplay,3000.,u_seri,u30)
+      call plevel(klon,klev,.false.,pplay,3000.,v_seri,v30)
+c 20 hPa
+      call plevel(klon,klev,.true. ,pplay,2000.,u_seri,u20)
+      call plevel(klon,klev,.false.,pplay,2000.,v_seri,v20)
+c 10 hPa
+      call plevel(klon,klev,.true. ,pplay,1000.,u_seri,u10)
+      call plevel(klon,klev,.false.,pplay,1000.,v_seri,v10)
+c
       call plevel(klon,klev,.true. ,pplay,50000.,zphi,phi500)
       call plevel(klon,klev,.true. ,paprs,50000.,omega,w500)
-
-cIM cf. FH     slp(:) = paprs(:,1)*exp(pphis(:)/(289.*t_seri(:,1)))
+c slp sea level pressure
       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
+      ENDDO
+c
+cIM sorties bilans energie cinetique et potentielle MJO
       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
-
+      DO i = 1, klon
+        d_u_oli(i,k) = d_u_oro(i,k) + d_u_lif(i,k)
+        d_v_oli(i,k) = d_v_oro(i,k) + d_v_lif(i,k)
+      ENDDO
+      ENDDO
 c=============================================================
 c   Ecriture des sorties
 c=============================================================
+#ifdef histREGDYN
+#include "write_histREGDYN.h"
+#endif
 
 #ifdef histISCCP
@@ -2739,4 +2484,9 @@
 #include "write_histday.h"
 #include "write_histmth.h"
+
+#ifdef histmthNMC
+#include "write_histmthNMC.h"
+#endif
+
 #include "write_histins.h"
 
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histhf.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histhf.h	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histhf.h	(revision 486)
@@ -54,4 +54,17 @@
       CALL histwrite(nid_hf,"phi500",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
 
+cIM cf FH
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, u_seri(:,1),zx_tmp_2d)
+      CALL histwrite(nid_hf,"u1",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, v_seri(:,1),zx_tmp_2d)
+      CALL histwrite(nid_hf,"v1",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
+      CALL histwrite(nid_hf,"cdrm",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragh,zx_tmp_2d)
+      CALL histwrite(nid_hf,"cdrh",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+
       if (ok_sync) then
         call histsync(nid_hf)
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h	(revision 485)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/write_histmth.h	(revision 486)
@@ -76,5 +76,7 @@
      .               ndex2d)
 c
-      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
+cIM: 071003
+      zx_tmp_fi2d(1:klon)=evap(1:klon)*86400.
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
       CALL histwrite(nid_mth,"evap",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
@@ -92,7 +94,12 @@
 c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollwdown,zx_tmp_2d)
-      CALL histwrite(nid_mth,"solldown",itau_w,zx_tmp_2d,iim*jjmp1,
+      CALL histwrite(nid_mth,"LWdnSFC",itau_w,zx_tmp_2d,iim*jjmp1,
      .               ndex2d)
-c
+cIM: 071003
+      zx_tmp_fi2d(1:klon)=sollw(1:klon)+sollwdown(1:klon)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d, zx_tmp_2d)
+      CALL histwrite(nid_mth,"LWupSFC",itau_w,zx_tmp_2d,iim*jjmp1,
+     .               ndex2d)
+cLWupSFC
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw0,zx_tmp_2d)
       CALL histwrite(nid_mth,"tops0",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
@@ -234,4 +241,20 @@
       CALL histwrite(nid_mth,"cldq",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
+      zx_tmp_fi2d(1:klon) = flwp(1:klon)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+      CALL histwrite(nid_mth,"lwp",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      zx_tmp_fi2d(1:klon) = fiwp(1:klon)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+      CALL histwrite(nid_mth,"iwp",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, flwc,zx_tmp_3d)
+      CALL histwrite(nid_mth,"lwcon",itau_w,zx_tmp_3d,
+     .               iim*jjmp1*klev,ndex3d)
+c
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, fiwc,zx_tmp_3d)
+      CALL histwrite(nid_mth,"iwcon",itau_w,zx_tmp_3d,
+     .               iim*jjmp1*klev,ndex3d)
+c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, ue,zx_tmp_2d)
       CALL histwrite(nid_mth,"ue",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
@@ -292,4 +315,10 @@
       CALL histwrite(nid_mth,"ovap",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
+cIM: 071003
+      zx_tmp_fi3d(1:klon,1:klev)=qx(1:klon,1:klev,ivap)/
+     .                         (1-qx(1:klon,1:klev,ivap))
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zx_tmp_fi3d, zx_tmp_3d)
+      CALL histwrite(nid_mth,"wvap",itau_w,zx_tmp_3d,
+     .                                   iim*jjmp1*klev,ndex3d)
 c
       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
@@ -333,4 +362,8 @@
      .                                   iim*jjmp1*klev,ndex3d)
 c
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t, zx_tmp_3d)
+      CALL histwrite(nid_mth,"dtphy",itau_w,zx_tmp_3d,
+     .                                   iim*jjmp1*klev,ndex3d)
+c
       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_dyn, zx_tmp_3d)
       CALL histwrite(nid_mth,"dtdyn",itau_w,zx_tmp_3d,
@@ -341,5 +374,7 @@
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_con, zx_tmp_3d)
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_con(1:klon,1:klev)/pdtphys
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtcon",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
@@ -349,6 +384,15 @@
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_lsc, zx_tmp_3d)
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_lsc(1:klon,1:klev)/pdtphys
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtlsc",itau_w,zx_tmp_3d,
+     .                                   iim*jjmp1*klev,ndex3d)
+cIM: 071003
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon, 1:klev)=(d_t_lsc(1:klon,1:klev)+
+     .                             d_t_eva(1:klon,1:klev))/pdtphys
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zx_tmp_fi3d, zx_tmp_3d)
+      CALL histwrite(nid_mth,"dtlschr",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
 c
@@ -357,5 +401,7 @@
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_vdf, zx_tmp_3d)
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_vdf(1:klon,1:klev)/pdtphys
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtvdf",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
@@ -365,5 +411,7 @@
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_eva, zx_tmp_3d)
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_eva(1:klon,1:klev)/pdtphys
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dteva",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
@@ -383,5 +431,7 @@
      .                                   iim*(jjmp1)*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_ajs, zx_tmp_3d)
+cIM: 101003 : K/30min ==> K/s
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtajs",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
@@ -391,21 +441,35 @@
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat, zx_tmp_3d)
+cIM: 101003 : K/day ==> K/s
+cIM: LMD_ARMIP3   zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)*pdtphys/RDAY 
+      zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY 
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtswr",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat0, zx_tmp_3d)
+cIM: 101003 : K/day ==> K/s      
+cIM: LMD_ARMIP3   zx_tmp_fi3d(1:klon,1:klev)=heat0(1:klon,1:klev)*pdtphys/RDAY
+      zx_tmp_fi3d(1:klon,1:klev)=heat0(1:klon,1:klev)/RDAY
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtsw0",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool, zx_tmp_3d)
+cIM: 101003 : K/day ==> K/s      
+cIM: LMD_ARMIP3     zx_tmp_fi3d(1:klon,1:klev)=-1.*cool(1:klon,1:klev)*pdtphys/RDAY
+      zx_tmp_fi3d(1:klon,1:klev)=-1.*cool(1:klon,1:klev)/RDAY
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtlwr",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool0, zx_tmp_3d)
+cIM: 101003 : K/day ==> K/s      
+cIM: LMD_ARMIP3     zx_tmp_fi3d(1:klon,1:klev)=-1.*cool0(1:klon,1:klev)*pdtphys/RDAY
+      zx_tmp_fi3d(1:klon,1:klev)=-1.*cool0(1:klon,1:klev)/RDAY
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtlw0",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
 c
-      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_ec, zx_tmp_3d)
+cIM: 101003 : deja en K/s      
+      zx_tmp_fi3d(1:klon,1:klev)=d_t_ec(1:klon,1:klev)
+      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
       CALL histwrite(nid_mth,"dtec",itau_w,zx_tmp_3d,
      .                                   iim*jjmp1*klev,ndex3d)
