--- trunk/phylmd/physiq.f	2014/07/02 18:39:15	99
+++ trunk/phylmd/physiq.f	2014/08/29 13:00:05	103
@@ -22,7 +22,7 @@
     USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &
          ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin
     USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &
-         ok_orodr, ok_orolf, soil_model
+         ok_orodr, ok_orolf
     USE clmain_m, ONLY: clmain
     use clouds_gno_m, only: clouds_gno
     USE comgeomphy, ONLY: airephy, cuphy, cvphy
@@ -242,7 +242,9 @@
     REAL, save:: fqsurf(klon, nbsrf)
     ! humidite de l'air au contact de la surface
 
-    REAL, save:: qsol(klon) ! hauteur d'eau dans le sol
+    REAL, save:: qsol(klon)
+    ! column-density of water in soil, in kg m-2
+
     REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
     REAL, save:: falbe(klon, nbsrf) ! albedo par type de surface
     REAL, save:: falblw(klon, nbsrf) ! albedo par type de surface
@@ -302,8 +304,11 @@
     REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
     REAL frac_nucl(klon, llm) ! idem (nucleation)
 
-    REAL, save:: rain_fall(klon) ! pluie
-    REAL, save:: snow_fall(klon) ! neige
+    REAL, save:: rain_fall(klon)
+    ! liquid water mass flux (kg/m2/s), positive down
+
+    REAL, save:: snow_fall(klon)
+    ! solid water mass flux (kg/m2/s), positive down
 
     REAL rain_tiedtke(klon), snow_tiedtke(klon)
 
@@ -391,7 +396,7 @@
     real zlongi
     REAL z_avant(klon), z_apres(klon), z_factor(klon)
     REAL za, zb
-    REAL zx_t, zx_qs, zdelta, zcor
+    REAL zx_t, zx_qs, zcor
     real zqsat(klon, llm)
     INTEGER i, k, iq, nsrf
     REAL, PARAMETER:: t_coup = 234.
@@ -738,7 +743,7 @@
     julien = MOD(NINT(rdayvrai), 360)
     if (julien == 0) julien = 360
 
-    forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg
+    forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
 
     ! Prescrire l'ozone :
     wo = ozonecm(REAL(julien), paprs)
@@ -792,15 +797,14 @@
     ! Couche limite:
 
     CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &
-         v_seri, julien, rmu0, co2_ppm, ftsol, soil_model, &
-         cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, &
-         fsnow, fqsurf, fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, &
-         fsolsw, fsollw, fder, rlat, frugs, firstcal, agesno, rugoro, &
-         d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, &
-         fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, &
-         u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, &
-         trmb3, plcl, fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, &
-         seaice)
+         v_seri, julien, rmu0, co2_ppm, ftsol, cdmmax, cdhmax, &
+         ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, &
+         fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, fsolsw, fsollw, &
+         fder, rlat, frugs, firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, &
+         d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &
+         q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, pblh, &
+         capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &
+         fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab)
 
     ! Incr\'ementation des flux
 
@@ -927,18 +931,13 @@
        dlw(i) = - 4. * RSIGMA * zxtsol(i)**3 
     ENDDO
 
-    ! Appeler la convection (au choix)
-
-    DO k = 1, llm
-       DO i = 1, klon
-          conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k) / dtphys
-          conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k) / dtphys
-       ENDDO
-    ENDDO
-
     IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
 
+    ! Appeler la convection (au choix)
+
     if (iflag_con == 2) then
+       conv_q = d_q_dyn + d_q_vdf / dtphys
+       conv_t = d_t_dyn + d_t_vdf / dtphys
        z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
        CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
             q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
@@ -963,26 +962,14 @@
        mfu = upwd + dnwd
        IF (.NOT. ok_gust) wd = 0.
 
-       ! Calcul des propri\'et\'es des nuages convectifs
-
-       DO k = 1, llm
-          DO i = 1, klon
-             IF (thermcep) THEN
-                zdelta = MAX(0., SIGN(1., rtt - t_seri(i, k)))
-                zqsat(i, k) = r2es * FOEEW(t_seri(i, k), zdelta) / play(i, k)
-                zqsat(i, k) = MIN(0.5, zqsat(i, k))
-                zqsat(i, k) = zqsat(i, k) / (1.-retv*zqsat(i, k))
-             ELSE
-                IF (t_seri(i, k) < t_coup) THEN
-                   zqsat(i, k) = qsats(t_seri(i, k))/play(i, k)
-                ELSE
-                   zqsat(i, k) = qsatl(t_seri(i, k))/play(i, k)
-                ENDIF
-             ENDIF
-          ENDDO
-       ENDDO
+       IF (thermcep) THEN
+          zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
+          zqsat = zqsat / (1. - retv * zqsat)
+       ELSE
+          zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
+       ENDIF
 
-       ! calcul des proprietes des nuages convectifs
+       ! Properties of convective clouds
        clwcon0 = fact_cldcon * clwcon0
        call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
             rnebcon0)
@@ -1221,8 +1208,7 @@
        DO i = 1, klon
           zx_t = t_seri(i, k)
           IF (thermcep) THEN
-             zdelta = MAX(0., SIGN(1., rtt-zx_t))
-             zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)
+             zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
              zx_qs = MIN(0.5, zx_qs)
              zcor = 1./(1.-retv*zx_qs)
              zx_qs = zx_qs*zcor