Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 383)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 385)
@@ -551,4 +551,5 @@
       save ok_newmicro
       save fact_cldcon,facttemps
+      real facteur
 
       integer iflag_cldcon
@@ -578,5 +579,5 @@
 c
       REAL t_seri(klon,klev), q_seri(klon,klev)
-      REAL ql_seri(klon,klev)
+      REAL ql_seri(klon,klev),qs_seri(klon,klev)
       REAL u_seri(klon,klev), v_seri(klon,klev)
 c
@@ -619,4 +620,21 @@
       character*40 vartitle
       character*20 varunits
+C     Variables liees au bilan d'energie et d'enthalpi
+      INTEGER   if_ebil ! level for energy conserv. dignostics
+      SAVE      if_ebil
+      REAL ztsol(klon)
+      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
+     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
+      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
+     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
+      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
+      REAL      d_h_vcol_phy
+      REAL      fs_bound, fq_bound
+      SAVE      d_h_vcol_phy
+      REAL      zero_v(klon)
+      CHARACTER*15 ztit
+      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
+      SAVE      ip_ebil
+      DATA      ip_ebil/2/
 c
 c Declaration des constantes et des fonctions thermodynamiques
@@ -627,4 +645,10 @@
 c======================================================================
       modname = 'physiq'
+      if_ebil = 2
+      IF (if_ebil.ge.1) THEN
+        DO i=1,klon
+          zero_v(i)=0.
+        END DO 
+      END IF 
       ok_sync=.TRUE.
       IF (nqmax .LT. 2) THEN
@@ -644,5 +668,6 @@
 c
        IF (debut) THEN
-
+C
+         IF (if_ebil.ge.1) d_h_vcol_phy=0.
 c
 c appel a la lecture du run.def physique
@@ -1526,4 +1551,12 @@
      .                "inst(X)", zsto,zout)
 c
+         CALL histdef(nid_ins, "cdrm", "Momentum drag coef.", "-",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "inst(X)", zsto,zout)
+c
+         CALL histdef(nid_ins, "cdrh", "Heat drag coef.", "-",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "inst(X)", zsto,zout)
+c
          CALL histdef(nid_ins, "rain", "Precipitation", "mm/day",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
@@ -1731,4 +1764,5 @@
          q_seri(i,k)  = qx(i,k,ivap)
          ql_seri(i,k) = qx(i,k,iliq)
+         qs_seri(i,k) = 0.
       ENDDO
       ENDDO
@@ -1748,5 +1782,27 @@
       ENDDO
       ENDIF
-c
+C
+      IF (if_ebil.ge.1) THEN 
+        DO i = 1, klon
+          ztsol(i) = 0.
+        ENDDO
+        DO i = 1, klon
+          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
+        ENDDO
+        ztit='after dynamic'
+        CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+C     Comme les tendances de la physique sont ajoute dans la dynamique,
+C     on devrait avoir que la variation d'entalpie par la dynamique
+C     est egale a la variation de la physique au pas de temps precedent.
+C     Donc la somme de ces 2 variations devrait etre nulle.
+        call diagphy(paire,ztit,ip_ebil
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
+     s      , fs_bound, fq_bound )
+      END IF 
+
 c Diagnostiquer la tendance dynamique
 c
@@ -1812,4 +1868,18 @@
       ENDDO
 c
+      IF (if_ebil.ge.2) THEN 
+        ztit='after reevap'
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(paire,ztit,ip_ebil
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+C
+      END IF 
+C
+c
 c Appeler la diffusion verticale (programme de couche limite)
 c
@@ -1825,5 +1895,6 @@
       DO nsrf = 1, nbsrf
       DO i = 1, klon
-         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
+c$$$         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001
+        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
       ENDDO
       ENDDO
@@ -1899,4 +1970,17 @@
       ENDDO
       ENDDO
+c
+      IF (if_ebil.ge.2) THEN 
+        ztit='after clmain'
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(paire,ztit,ip_ebil
+     e      , zero_v, zero_v, zero_v, zero_v, sens
+     e      , evap  , zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+C
 c
 c Incrementer la temperature du sol
@@ -2072,4 +2156,17 @@
         ENDDO
       ENDDO
+c
+      IF (if_ebil.ge.2) THEN 
+        ztit='after convect'
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(paire,ztit,ip_ebil
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, rain_con, snow_con, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+C
       IF (check) THEN
           za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
@@ -2130,4 +2227,11 @@
       ENDDO
       ENDDO
+c
+      IF (if_ebil.ge.2) THEN 
+        ztit='after dry_adjust'
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+      END IF 
 
 
@@ -2168,6 +2272,7 @@
 c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
 c   relaxation des ratqs
-         facttemps=exp(-pdtphys/1.e4)
-         ratqs(:,:)=max(ratqs(:,:)*facttemps,ratqss(:,:))
+c         facttemps=exp(-pdtphys/1.e4)
+         facteur=exp(-pdtphys*facttemps)
+         ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
          ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
 c         print*,'calcul des ratqs fini'
@@ -2214,4 +2319,16 @@
       ENDIF
 c
+      IF (if_ebil.ge.2) THEN 
+        ztit='after fisrt'
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(paire,ztit,ip_ebil
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, rain_lsc, snow_lsc, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+c
 c-------------------------------------------------------------------
 c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
@@ -2239,8 +2356,9 @@
 c  convection et du calcul du pas de temps précédent diminué d'un facteur
 c  facttemps
-      facttemps=pdtphys/1.e4
+c      facttemps=pdtphys/1.e4
+      facteur = pdtphys *facttemps
       do k=1,klev
          do i=1,klon
-            rnebcon(i,k)=rnebcon(i,k)*facttemps
+            rnebcon(i,k)=rnebcon(i,k)*facteur
             if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
      s      then
@@ -2278,4 +2396,11 @@
          snow_fall(i) = snow_con(i) + snow_lsc(i)
       ENDDO
+c
+      IF (if_ebil.ge.2) THEN 
+        ztit="after diagcld"
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+      END IF 
 c
 c Calculer l'humidite relative pour diagnostique
@@ -2356,4 +2481,17 @@
       ENDDO
 c
+      IF (if_ebil.ge.2) THEN 
+        ztit='after rad'
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(paire,ztit,ip_ebil
+     e      , topsw, toplw, solsw, sollw, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+c
+c
 c Calculer l'hydrologie de la surface
 c
@@ -2457,4 +2595,12 @@
 c
       ENDIF ! fin de test sur ok_orolf
+c
+      IF (if_ebil.ge.2) THEN 
+        ztit='after orography'
+        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+      END IF 
+c
 c
 cAA
@@ -2520,5 +2666,23 @@
 c
 c
-
+      IF (if_ebil.ge.1) THEN 
+        ztit='after physic'
+        CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+C     Comme les tendances de la physique sont ajoute dans la dynamique,
+C     on devrait avoir que la variation d'entalpie par la dynamique
+C     est egale a la variation de la physique au pas de temps precedent.
+C     Donc la somme de ces 2 variations devrait etre nulle.
+        call diagphy(paire,ztit,ip_ebil
+     e      , topsw, toplw, solsw, sollw, sens
+     e      , evap, rain_fall, snow_fall, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+C
+      d_h_vcol_phy=d_h_vcol
+C
+      END IF 
+C
       IF (ok_journe) THEN
 c
@@ -3272,4 +3436,9 @@
       CALL histwrite(nid_ins,"snow",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
 
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
+      CALL histwrite(nid_ins,"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_ins,"cdrh",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
