--- trunk/Sources/phylmd/clmain.f	2017/11/10 15:16:48	239
+++ trunk/Sources/phylmd/clmain.f	2017/11/14 14:38:36	243
@@ -106,7 +106,7 @@
     ! dflux_q derive du flux latent
     ! IM "slab" ocean
 
-    REAL, intent(out):: ycoefh(:, :) ! (klon, klev)
+    REAL, intent(out):: ycoefh(:, 2:) ! (klon, 2:klev)
     ! Pour pouvoir extraire les coefficients d'\'echange, le champ
     ! "ycoefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
     ! ce champ sur les quatre sous-surfaces du mod\`ele.
@@ -169,7 +169,7 @@
     REAL yu(klon, klev), yv(klon, klev)
     REAL yt(klon, klev), yq(klon, klev)
     REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
-    REAL ycoefm0(klon, klev), ycoefh0(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)
@@ -310,18 +310,15 @@
              END DO
           END DO
 
-          ! calculer Cdrag et les coefficients d'echange
           CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &
                yrugos, yu, yv, yt, yq, yqsurf(:knon), coefm(:knon, :), &
                coefh(:knon, :), ycdragm(:knon), ycdragh(:knon))
 
           IF (iflag_pbl == 1) THEN
-             CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, 2:), &
-                  ycoefh0(:knon, 2:))
-             ycoefm0(:knon, 1) = 0.
-             ycoefh0(:knon, 1) = 0.
-             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, 2:))
-             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))
+             CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0(:knon, :), &
+                  ycoefh0(:knon, :))
+             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
+             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
              ycdragm(:knon) = max(ycdragm(:knon), 0.)
              ycdragh(:knon) = max(ycdragh(:knon), 0.)
           END IF
@@ -335,10 +332,10 @@
           IF (ok_kzmin) THEN
              ! Calcul d'une diffusion minimale pour les conditions tres stables
              CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
-                  ycdragm(:knon), ycoefh0(:knon, 2:))
-             ycoefm0(:knon, 2:) = ycoefh0(:knon, 2:)
-             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, 2:))
-             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, 2:))
+                  ycdragm(:knon), ycoefh0(:knon, :))
+             ycoefm0(:knon, :) = ycoefh0(:knon, :)
+             coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
+             coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
           END IF
 
           IF (iflag_pbl >= 6) THEN
@@ -418,18 +415,6 @@
              y_dflux_q(j) = y_dflux_q(j) * ypct(j)
           END DO
 
-          DO k = 2, klev
-             DO j = 1, knon
-                i = ni(j)
-                coefh(j, k) = coefh(j, k) * ypct(j)
-                coefm(j, k) = coefm(j, k) * ypct(j)
-             END DO
-          END DO
-          DO j = 1, knon
-             i = ni(j)
-             ycdragh(j) = ycdragh(j) * ypct(j)
-             ycdragm(j) = ycdragm(j) * ypct(j)
-          END DO
           DO k = 1, klev
              DO j = 1, knon
                 i = ni(j)
@@ -466,8 +451,8 @@
              agesno(i, nsrf) = yagesno(j)
              fqcalving(i, nsrf) = y_fqcalving(j)
              ffonte(i, nsrf) = y_ffonte(j)
-             cdragh(i) = cdragh(i) + ycdragh(j)
-             cdragm(i) = cdragm(i) + ycdragm(j)
+             cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
+             cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
              dflux_t(i) = dflux_t(i) + y_dflux_t(j)
              dflux_q(i) = dflux_q(i) + y_dflux_q(j)
           END DO
@@ -492,18 +477,9 @@
                 d_v(i, k) = d_v(i, k) + y_d_v(j, k)
              END DO
           END DO
-          
-          DO j = 1, knon
-             i = ni(j)
-             DO k = 2, klev
-                ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
-             END DO
-          END DO
 
-          DO j = 1, knon
-             i = ni(j)
-             ycoefh(i, 1) = ycoefh(i, 1) + ycdragh(j)
-          END DO
+          forall (k = 2:klev) ycoefh(ni(:knon), k) &
+               = ycoefh(ni(:knon), k) + coefh(:knon, k) * ypct(:knon)
 
           ! diagnostic t, q a 2m et u, v a 10m