Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 392)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 393)
@@ -482,5 +482,4 @@
       DO j = 1, knon
          yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
-     $      + 0.000014 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))
          yrugm(j) = MAX(1.5e-05,yrugm(j))
       ENDDO
@@ -710,4 +709,6 @@
       real emis_new(klon), z0_new(klon)
       real pctsrf_new(klon,nbsrf)
+c JLD
+      real zzpk
       
 c
@@ -773,6 +774,7 @@
          zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
 c
-         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
-         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
+         zzpk=(pplay(i,klev)/psref(i))**RKAPPA
+         zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)
+         zx_ch(i,klev) = (local_h(i,klev)*zzpk*delp(i,klev)
      .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
          zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
@@ -788,7 +790,8 @@
          zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
 c
-         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
+         zzpk=(pplay(i,k)/psref(i))**RKAPPA
+         zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k)
      .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
-         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
+         zx_ch(i,k) = (local_h(i,k)*zzpk*delp(i,k)
      .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
      .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
@@ -810,6 +813,7 @@
          zx_dq(i,1) = -1. * RG / zx_buf1(i)
 c
-         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
-         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
+         zzpk=(pplay(i,1)/psref(i))**RKAPPA
+         zx_buf2(i) = zzpk*delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
+         zx_ch(i,1) = (local_h(i,1)*zzpk*delp(i,1)
      .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
      .                /zx_buf2(i)
@@ -1101,7 +1105,5 @@
       REAL zcfm1(klon), zcfm2(klon)
       REAL zcfh1(klon), zcfh2(klon)
-c$$$      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
-c$$$cPB differenciation coefficient de frottement drag et flux chaleur
-      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn,zcdh, rugh
+      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
       REAL zscf, zucf, zcr
       REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
@@ -1186,126 +1188,6 @@
 c Calculer le frottement au sol (Cdrag)
 c
-c$$$c PB
-c$$$c essais d'itération pour l'océan
-c
-      rugh = 1.3 e-4
-      IF (nsrf.EQ.is_oce) THEN 
-      DO k=1,10
-c$$$        WRITE(*,*) 'k',k
-c$$$        WRITE(*,*) rugos(100)
-      DO i = 1, knon
-         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
-         zdphi=zgeop(i,1)
-         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
-c$$$         ztsolv = ts(i) * (1.0+RETV*qsurf(i))
-         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
-     .       *(1.+RETV*q(i,1))
-         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
-         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
-c PB ajout drag neutre pour flux chaleur
-         zcdh = ckap**2/(log(1.+zgeop(i,1)/(RG*rugos(i)))
-     $       * log(1.+zgeop(i,1)/(RG*rugh)))
-         IF (zri(i) .ge. 0.) THEN ! situation stable
-           IF (.NOT.zxli) THEN
-           zscf=SQRT(1.+cd*ABS(zri(i)))
-           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
-!           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf)
-           zcfm1(i) = zcdn * FRIV
-           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
-!           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
-c PB avec drag neutre pour flux chaleur different
-c           zcfh1(i) = zcdn * FRIH
-           zcfh1(i) = zcdh * FRIH
-           pcfm(i,1) = zcfm1(i)
-           pcfh(i,1) = zcfh1(i)
-           ELSE
-           pcfm(i,1) = zcdn* fsta(zri(i))
-           pcfh(i,1) = zcdn* fsta(zri(i))
-           ENDIF
-         ELSE ! situation instable
-           IF (.NOT.zxli) THEN
-           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
-     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
-           zcfm2(i) = zcdn*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
-c PB ajout pour prendre drag neutre des flux chaleur different drag neutre vent
-           zucf=1./(1.+3.0*cb*cc*zcdh*SQRT(ABS(zri(i))
-     .              *(1.0+zgeop(i,1)/(RG*rugh))))
-c$$$           zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
-           zcfh2(i) = zcdh*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
-           pcfm(i,1) = zcfm2(i)
-           pcfh(i,1) = zcfh2(i)
-           ELSE
-           pcfm(i,1) = zcdn* fins(zri(i))
-           pcfh(i,1) = zcdn* fins(zri(i))
-           ENDIF
-           zcr = (0.0016/(zcdh*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
-           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdh*(1.0+zcr**1.25)**(1./1.25)
-         ENDIF
-C
-c$$$C PB test drag
-c$$$         pcfm(i,1)=zcdn
-c$$$         pcfh(i,1) = zcdn
-         rugos(i)= 0.018*pcfm(i,1) * zdu2/RG 
-     $       +0.11*0.000014/sqrt(pcfm(i,1) * zdu2)
-         rugh = 0.62*0.000014/sqrt(pcfm(i,1) * zdu2)+1.4e-5
-      END DO
-      END DO 
-
-      ELSE 
-      DO i = 1, knon
-         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
-         zdphi=zgeop(i,1)
-         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
-c         ztsolv = ts(i) * (1.0+RETV*qsurf(i))
-         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
-     .       *(1.+RETV*q(i,1))
-         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
-         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
-c pb ajout pour avoir drage neutre flux chaleur differents
-         zcdh = ckap**2/(log(1.+zgeop(i,1)/(RG*rugos(i)))
-     $       * log(1.+zgeop(i,1)/(RG*1.3e-4)))
-         IF (zri(i) .ge. 0.) THEN ! situation stable
-           IF (.NOT.zxli) THEN
-           zscf=SQRT(1.+cd*ABS(zri(i)))
-           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
-!           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf)
-           zcfm1(i) = zcdn * FRIV
-           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
-!           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
-c$$$C Modif pour drag neutre flux chaleur differents
-c$$$           zcfh1(i) = zcdn * FRIH
-           zcfh1(i) = zcdh * FRIH
-           pcfm(i,1) = zcfm1(i)
-           pcfh(i,1) = zcfh1(i)
-           ELSE
-           pcfm(i,1) = zcdn* fsta(zri(i))
-           pcfh(i,1) = zcdn* fsta(zri(i))
-           ENDIF
-         ELSE ! situation instable
-           IF (.NOT.zxli) THEN
-           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
-     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
-           zcfm2(i) = zcdn*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
-C PB ajout pour drage neutre flux chaleur differents
-           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
-     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
-c$$$           zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
-           zcfh2(i) = zcdh*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
-           pcfm(i,1) = zcfm2(i)
-           pcfh(i,1) = zcfh2(i)
-           ELSE
-           pcfm(i,1) = zcdn* fins(zri(i))
-           pcfh(i,1) = zcdn* fins(zri(i))
-           ENDIF
-           zcr = (0.0016/(zcdn*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
-           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
-         ENDIF
-C
-c$$$C PB test drag
-c$$$         pcfm(i,1)=zcdn
-c$$$         pcfh(i,1) = zcdn
-      END DO
-c$$$CPB fin test iterations
-      ENDIF 
+      CALL clcdrag(knon, nsrf, zxli, u, v, t, q, zgeop,
+     .             ts, qsurf, rugos, pcfm, pcfh, zcdn, zri)
 c
 c Calculer les coefficients turbulents dans l'atmosphere
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F	(revision 392)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F	(revision 393)
@@ -104,4 +104,7 @@
       REAL zprec_cond(klon)
 cAA
+      REAL zmair, zcpair, zcpeau
+C     Pour la conversion eau-neige
+      REAL zlh_solid(klon), zm_solid
 c---------------------------------------------------------------
 c
@@ -209,6 +212,25 @@
       ENDDO
 c
+c Calculer la varition de temp. de l'air du a la chaleur sensible
+C transporter par la pluie.
+C Il resterait a rajouter cet effet de la chaleur sensible sur les
+C flux de surface, du a la diff. de temp. entre le 1er niveau et la
+C surface.
+C
+      DO i = 1, klon
+        zmair=(paprs(i,k)-paprs(i,k+1))/RG
+        zcpair=RCPD*(1.0+RVTMP2*zq(i))
+        zcpeau=RCPD*RVTMP2
+        zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau
+     $      + zmair*zcpair*zt(i) )
+     $      / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
+CC        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
+      ENDDO
+c
+c
 c Calculer l'evaporation de la precipitation
 c
+
+
       IF (evap_prec) THEN
       DO i = 1, klon
@@ -363,5 +385,6 @@
       DO i = 1, klon
          zq(i) = zq(i) - zcond(i)
-         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
+c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
+         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
       ENDDO
 c
@@ -484,8 +507,21 @@
       IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
          snow(i) = zrfl(i)
+         zlh_solid(i) = RLSTT-RLVTT
       ELSE
          rain(i) = zrfl(i)
-      ENDIF
-      ENDDO
+         zlh_solid(i) = 0.
+      ENDIF
+      ENDDO
+C
+C For energy conservation : when snow is present, the solification
+c latent heat is considered.
+      DO k = 1, klev
+        DO i = 1, klon
+          zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
+          zmair=(paprs(i,k)-paprs(i,k+1))/RG
+          zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
+          d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
+        END DO 
+      END DO
 c
       RETURN
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 392)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 393)
@@ -39,4 +39,5 @@
 #include "dimphy.h"
 #include "raddim.h"
+#include "YOETHF.h"
 c
       real rmu0(klon), fract(klon), dist
@@ -88,5 +89,5 @@
       REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
       REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
-
+      REAL*8 zznormcp
 c
 c-------------------------------------------
@@ -215,9 +216,18 @@
       ENDDO
       DO k = 1, kflev
+c      DO i = 1, kdlon
+c         heat(iof+i,k) = zheat(i,k)
+c         cool(iof+i,k) = zcool(i,k)
+c         heat0(iof+i,k) = zheat0(i,k)
+c         cool0(iof+i,k) = zcool0(i,k)
+c      ENDDO
       DO i = 1, kdlon
-         heat(iof+i,k) = zheat(i,k)
-         cool(iof+i,k) = zcool(i,k)
-         heat0(iof+i,k) = zheat0(i,k)
-         cool0(iof+i,k) = zcool0(i,k)
+C        scale factor to take into account the difference between
+C        dry air and watter vapour scpecific heat capacity
+         zznormcp=1.0+RVTMP2*PWV(i,k)
+         heat(iof+i,k) = zheat(i,k)/zznormcp
+         cool(iof+i,k) = zcool(i,k)/zznormcp
+         heat0(iof+i,k) = zheat0(i,k)/zznormcp
+         cool0(iof+i,k) = zcool0(i,k)/zznormcp
       ENDDO
       ENDDO
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/suphec.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/suphec.F	(revision 392)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/suphec.F	(revision 393)
@@ -123,5 +123,5 @@
 C              ---------------------------------------------
 C
-      RCW=4218.
+      RCW=RCPV
       WRITE(UNIT=6,FMT='('' *** Thermodynamic, liquid  ***'')')
       WRITE(UNIT=6,FMT='(''         Cw   = '',E13.7)') RCW
@@ -132,5 +132,5 @@
 C              --------------------------------------------
 C
-      RCS=2106.
+      RCS=RCPV
       WRITE(UNIT=6,FMT='('' *** thermodynamic, solid   ***'')')
       WRITE(UNIT=6,FMT='(''         Cs   = '',E13.7)') RCS
