--- trunk/Sources/phylmd/clmain.f	2017/11/14 14:56:42	244
+++ trunk/Sources/phylmd/clmain.f	2018/01/05 16:40:13	248
@@ -21,6 +21,7 @@
     ! ne tient pas compte de la diff\'erentiation des sous-fractions
     ! de sol.
 
+    use clcdrag_m, only: clcdrag
     use clqh_m, only: clqh
     use clvent_m, only: clvent
     use coefkz_m, only: coefkz
@@ -150,7 +151,8 @@
     real y_run_off_lic_0(klon)
     REAL rugmer(klon)
     REAL ytsoil(klon, nsoilmx)
-    REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
+    REAL yts(klon), ypct(klon), yz0_new(klon)
+    real yrugos(klon) ! longeur de rugosite (en m)
     REAL yalb(klon)
     REAL snow(klon), yqsurf(klon), yagesno(klon)
     real yqsol(klon) ! column-density of water in soil, in kg m-2
@@ -171,7 +173,6 @@
     REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
     REAL ycoefm0(klon, 2:klev), ycoefh0(klon, 2:klev)
     REAL yzlay(klon, klev), zlev(klon, klev + 1), yteta(klon, klev)
-    REAL ykmm(klon, klev + 1), ykmn(klon, klev + 1)
     REAL yq2(klon, klev + 1)
     REAL delp(klon, klev)
     INTEGER i, k, nsrf
@@ -201,6 +202,7 @@
 
     REAL qairsol(klon), zgeo1(klon)
     REAL rugo1(klon)
+    REAL zgeop(klon, klev)
 
     !------------------------------------------------------------
 
@@ -310,9 +312,25 @@
              END DO
           END DO
 
-          CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &
-               yrugos, yu, yv, yt, yq, yqsurf(:knon), ycoefm(:knon, :), &
-               ycoefh(:knon, :), ycdragm(:knon), ycdragh(:knon))
+          ! Calculer les géopotentiels de chaque couche:
+
+          zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
+               + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
+
+          DO k = 2, klev
+             zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
+                  * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
+                  * (ypplay(:knon, k - 1) - ypplay(:knon, k))
+          ENDDO
+
+          CALL clcdrag(nsrf, yu(:knon, 1), yv(:knon, 1), yt(:knon, 1), &
+               yq(:knon, 1), zgeop(:knon, 1), yts(:knon), yqsurf(:knon), &
+               yrugos(:knon), ycdragm(:knon), ycdragh(:knon)) 
+
+          CALL coefkz(nsrf, ypaprs(:knon, :), ypplay(:knon, :), ksta, &
+               ksta_ter, yts(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
+               yq(:knon, :), zgeop(:knon, :), ycoefm(:knon, :), &
+               ycoefh(:knon, :))
 
           IF (iflag_pbl == 1) THEN
              CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, :), &
@@ -375,9 +393,7 @@
              ustar(:knon) = ustarhb(yu(:knon, 1), yv(:knon, 1), ycdragm(:knon))
              CALL yamada4(dtime, rg, zlev(:knon, :), yzlay(:knon, :), &
                   yu(:knon, :), yv(:knon, :), yteta(:knon, :), yq2(:knon, :), &
-                  ykmm(:knon, :), ykmn(:knon, :), ustar(:knon))
-             ycoefm(:knon, :) = ykmm(:knon, 2:klev)
-             ycoefh(:knon, :) = ykmn(:knon, 2:klev)
+                  ycoefm(:knon, :), ycoefh(:knon, :), ustar(:knon))
           END IF
 
           CALL clvent(dtime, yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
@@ -504,7 +520,7 @@
 
           CALL stdlevvar(klon, knon, nsrf, u1(:knon), v1(:knon), tair1(:knon), &
                qair1, zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, &
-               yq2m, yt10m, yq10m, wind10m(:knon), ustar)
+               yq2m, yt10m, yq10m, wind10m(:knon), ustar(:knon))
 
           DO j = 1, knon
              i = ni(j)