Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 86)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 98)
@@ -1,6 +1,8 @@
       SUBROUTINE clmain(dtime,pctsrf,t,q,u,v,
-     .                  soil_model,ts,soilcap,soilflux,
-     .                  paprs,pplay,radsol,snow,qsol,
-     .                  xlat, rugos,
+     .                  ok_veget,ts,
+     .                  paprs,pplay,radsol,snow,qsol,evap,albe,
+     .                  rain_f, snow_f, solsw, sollw,
+     .                  rlon, rlat, rugos,
+     .                  debut, lafin,
      .                  d_t,d_q,d_u,d_v,d_ts,
      .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
@@ -41,5 +43,5 @@
 c beta-----input-R- coefficient de l'evaporation reelle (0 a 1)
 c dif_grnd-input-R- coeff. de diffusion (chaleur) vers le sol profond
-c xlat-----input-R- latitude en degree
+c rlat-----input-R- latitude en degree
 c rugos----input-R- longeur de rugosite (en m)
 c
@@ -75,12 +77,13 @@
       REAL u(klon,klev), v(klon,klev)
       REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
-      REAL xlat(klon)
+      REAL rlon(klon), rlat(klon)
       REAL d_t(klon, klev), d_q(klon, klev)
       REAL d_u(klon, klev), d_v(klon, klev)
-      REAL flux_t(klon,klev), flux_q(klon,klev)
+      REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)
       REAL dflux_t(klon), dflux_q(klon)
-      REAL flux_u(klon,klev), flux_v(klon,klev)
+      REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)
       REAL rugmer(klon)
       REAL cdragh(klon), cdragm(klon)
+      LOGICAL debut, lafin, ok_veget
 cAA      INTEGER itr
 cAA      REAL tr(klon,klev,nbtr)
@@ -93,4 +96,8 @@
       REAL snow(klon,nbsrf)
       REAL qsol(klon,nbsrf)
+      REAL evap(klon,nbsrf)
+      REAL albe(klon,nbsrf)
+      real rain_f(klon), snow_f(klon)
+      REAL sollw(klon), solsw(klon)
       REAL rugos(klon,nbsrf)
 cAA
@@ -103,6 +110,9 @@
 c======================================================================
       REAL yts(klon), yrugos(klon), ypct(klon)
-      REAL ycal(klon), ybeta(klon), ydif(klon)
+      REAL ycal(klon), ybeta(klon), ydif(klon), yalb(klon),yevap(klon)
       REAL yu1(klon), yv1(klon)
+      real ysnow(klon), yqsol(klon)
+      real yrain_f(klon), ysnow_f(klon)
+      real ysollw(klon), ysolsw(klon), ysolswnet(klon)
       REAL yrugm(klon), yrads(klon)
       REAL y_d_ts(klon)
@@ -144,6 +154,4 @@
 c======================================================================
 
-      write(*,*)'CLMAIN.NEW'
-
       DO k = 1, klev   ! epaisseur de couche
       DO i = 1, klon
@@ -174,15 +182,20 @@
          d_ts(i,nsrf) = 0.0
       ENDDO
-      ENDDO
+      END DO
+C§§§ PB
+      flux_t = 0.
+      flux_q = 0.
+      flux_u = 0.
+      flux_v = 0.
       DO k = 1, klev
       DO i = 1, klon
          d_t(i,k) = 0.0
          d_q(i,k) = 0.0
-         flux_t(i,k) = 0.0
-         flux_q(i,k) = 0.0
+c$$$         flux_t(i,k) = 0.0
+c$$$         flux_q(i,k) = 0.0
          d_u(i,k) = 0.0
          d_v(i,k) = 0.0
-         flux_u(i,k) = 0.0
-         flux_v(i,k) = 0.0
+c$$$         flux_u(i,k) = 0.0
+c$$$         flux_v(i,k) = 0.0
          zcoefh(i,k) = 0.0
       ENDDO
@@ -203,21 +216,23 @@
 c
 c prescrire les proprietes du sol:
-      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
-      IF (.NOT.soil_model) THEN
-         DO i = 1, klon
-            cal(i) = RCPD * capsol(i)
-            totalflu(i) = radsol(i)
-         ENDDO
-      ELSE
-         DO i = 1, klon
-            totalflu(i) = soilflux(i,nsrf) + radsol(i)
-            IF (nsrf.EQ.is_oce) THEN
-               cal(i) = 0.0
-            ELSE
-               cal(i) = RCPD / soilcap(i,nsrf)
-            ENDIF
-         ENDDO
-      ENDIF
-c
+c      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
+c      IF (.NOT.soil_model) THEN
+c         DO i = 1, klon
+c            cal(i) = RCPD * capsol(i)
+c            totalflu(i) = radsol(i)
+c         ENDDO
+c      ELSE
+c         DO i = 1, klon
+c            totalflu(i) = soilflux(i,nsrf) + radsol(i)
+c            IF (nsrf.EQ.is_oce) THEN
+c               cal(i) = 0.0
+c            ELSE
+c               cal(i) = RCPD / soilcap(i,nsrf)
+c            ENDIF
+c         ENDDO
+c      ENDIF
+c
+      totalflu = radsol
+
 c chercher les indices:
       DO j = 1, klon
@@ -237,11 +252,19 @@
         ypct(j) = pctsrf(i,nsrf)
         yts(j) = ts(i,nsrf)
+        ysnow(j) = snow(i,nsrf)
+        yevap(j) = evap(i,nsrf)
+        yqsol(j) = qsol(i,nsrf)
+        yalb(j) = albe(i,nsrf)
+        yrain_f(j) = rain_f(i)
+        ysnow_f(j) = snow_f(i)
+        ysolsw(j) = solsw(i)
+        ysollw(j) = sollw(i)
         yrugos(j) = rugos(i,nsrf)
         yu1(j) = u1lay(i)
         yv1(j) = v1lay(i)
         yrads(j) = totalflu(i)
-        ycal(j) = cal(i)
-        ybeta(j) = beta(i)
-        ydif(j) = dif_grnd(i)
+c        ycal(j) = cal(i)
+c        ybeta(j) = beta(i)
+c        ydif(j) = dif_grnd(i)
         ypaprs(j,klev+1) = paprs(i,klev+1)
       ENDDO
@@ -259,18 +282,4 @@
       ENDDO
 c
-cAA      IF (itr.GE.1) THEN
-cAA      DO it = 1, itr
-cAA        DO k = 1, klev
-cAA        DO j = 1, knon
-cAA        i = ni(j)
-cAA          ytr(j,k,it) = tr(i,k,it)
-cAA        ENDDO
-cAA        ENDDO
-cAA        DO j = 1, knon
-cAA        i = ni(j)
-cAA          yflxsrf(j,it) = flux_surf(i,it)
-cAA        ENDDO
-cAA      ENDDO
-cAA      ENDIF
 c
 c calculer Cdrag et les coefficients d'echange
@@ -287,45 +296,11 @@
       ENDDO
 c
-c parametrisation non-locale:
-      IF (ok_nonloc) THEN
-         DO i = 1, knon
-            y_cd_h(i) = ycoefh(i,1)
-            y_cd_m(i) = ycoefm(i,1)
-         ENDDO
-         CALL nonlocal(knon, ypaprs, ypplay,
-     .        yts,ybeta,yu,yv,yt,yq,
-     .        y_cd_h, y_cd_m, ycoefm0, ycoefh0, ygamt, ygamq)
-         DO k = 1, klev
-         DO i = 1, knon
-            ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
-            ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
-         ENDDO
-         ENDDO
-      ELSE
-         IF (.NOT.contreg) THEN
-            DO k = 2, klev
-            DO i = 1, knon
-               ygamq(i,k) = 0.0
-               ygamt(i,k) = 0.0
-            ENDDO
-            ENDDO
-         ELSE
-         DO k = 3, klev
-         DO i = 1, knon
-            ygamq(i,k) = 0.0
-            ygamt(i,k) = -1.0E-03
-         ENDDO
-         ENDDO
-         DO i = 1, knon
-            ygamq(i,2) = 0.0
-            ygamt(i,2) = -2.5E-03
-         ENDDO
-         ENDIF
-      ENDIF
 c
 c calculer la diffusion de "q" et de "h"
-      CALL clqh(knon, dtime, nsrf,yu1, yv1,
+      CALL clqh(knon, dtime, nsrf, ni, pctsrf, rlon, rlat,
+     e          yu1, yv1,
      e          ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
-     e          ycal,ybeta,ydif,ygamt,ygamq,
+     e          yevap,yalb, ysnow, yqsol, yrain_f, ysnow_f,
+     e          ysollw, ysolsw,
      s          y_d_t, y_d_q, y_d_ts, 
      s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
@@ -363,21 +338,32 @@
 c
       DO k = 1, klev
+        DO j = 1, knon
+          i = ni(j)
+          ycoefh(j,k) = ycoefh(j,k) * ypct(j)
+          ycoefm(j,k) = ycoefm(j,k) * ypct(j)
+          y_d_t(j,k) = y_d_t(j,k) * ypct(j)
+          y_d_q(j,k) = y_d_q(j,k) * ypct(j)
+C§§§ PB
+          flux_t(i,k,nsrf) = y_flux_t(j,k)
+          flux_q(i,k,nsrf) = y_flux_q(j,k)
+          flux_u(i,k,nsrf) = y_flux_u(j,k)
+          flux_v(i,k,nsrf) = y_flux_v(j,k)
+c$$$ PB        y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
+c$$$ PB        y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
+          y_d_u(j,k) = y_d_u(j,k) * ypct(j)
+          y_d_v(j,k) = y_d_v(j,k) * ypct(j)
+c$$$ PB        y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
+c$$$ PB        y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
+        ENDDO
+      ENDDO
+
+      evap(:,nsrf) = - flux_q(:,1,nsrf)
+c
       DO j = 1, knon
-         ycoefh(j,k) = ycoefh(j,k) * ypct(j)
-         ycoefm(j,k) = ycoefm(j,k) * ypct(j)
-         y_d_t(j,k) = y_d_t(j,k) * ypct(j)
-         y_d_q(j,k) = y_d_q(j,k) * ypct(j)
-         y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
-         y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
-         y_d_u(j,k) = y_d_u(j,k) * ypct(j)
-         y_d_v(j,k) = y_d_v(j,k) * ypct(j)
-         y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
-         y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
-      ENDDO
-      ENDDO
-c
-      DO j = 1, knon
-      i = ni(j)
+         i = ni(j)
          d_ts(i,nsrf) = y_d_ts(j)
+         albe(i,nsrf) = yalb(j)
+         snow(i,nsrf) = ysnow(j)
+         qsol(i,nsrf) = yqsol(j)
          rugmer(i) = yrugm(j)
          cdragh(i) = cdragh(i) + ycoefh(j,1)
@@ -400,10 +386,10 @@
          d_t(i,k) = d_t(i,k) + y_d_t(j,k)
          d_q(i,k) = d_q(i,k) + y_d_q(j,k)
-         flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
-         flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
+c$$$ PB        flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
+c$$$         flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
          d_u(i,k) = d_u(i,k) + y_d_u(j,k)
          d_v(i,k) = d_v(i,k) + y_d_v(j,k)
-         flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
-         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
+c$$$  PB       flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
+c$$$         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
          zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
       ENDDO
@@ -430,9 +416,13 @@
 99999 CONTINUE
 c
+
       RETURN
       END
-      SUBROUTINE clqh(knon,dtime,nisurf,u1lay,v1lay,coef,
+      SUBROUTINE clqh(knon,dtime,nisurf,knindex,pctsrf, rlon, rlat,
+     e                u1lay,v1lay,coef,
      e                t,q,ts,paprs,pplay,
-     e                delp,radsol,cal,beta,dif_grnd, gamt,gamq,
+     e                delp,radsol,evap,albedo,snow,qsol, 
+     e                precip_rain, precip_snow,
+     e                lwdown, swdown,
      s                d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
 
@@ -446,4 +436,8 @@
 #include "dimensions.h"
 #include "dimphy.h"
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+#include "indicesol.h"
 c Arguments:
       INTEGER knon
@@ -462,8 +456,16 @@
       REAL q(klon,klev)       ! humidite specifique (kg/kg)
       REAL ts(klon)           ! temperature du sol (K)
+      REAL evap(klon)         ! evaporation au sol
       REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
       REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
       REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
       REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
+      REAL albedo(klon)       ! albedo de la surface
+      REAL snow(klon)         ! hauteur de neige
+      REAL qsol(klon)         ! humidite de la surface
+      real precip_rain(klon), precip_snow(klon)
+      integer knindex(klon)
+      real pctsrf(klon,nbsrf)
+      real rlon(klon), rlat(klon)
 c
       REAL d_t(klon,klev)     ! incrementation de "t"
@@ -521,8 +523,6 @@
 c======================================================================
       REAL zcor, zdelta, zcvm5
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "indicesol.h"
+      logical contreg
+      parameter (contreg=.true.)
 c======================================================================
 c Rajout pour l'interface
@@ -530,13 +530,10 @@
       integer jour
       integer nisurf
-      integer knindex(klon)
       logical debut, lafin, ok_veget
-      real rlon(klon), rlat(klon)
-      real zlev(klon), zlflu(klon)
+      real zlev1(klon)
       real temp_air(klon), spechum(klon)
       real hum_air(klon), ccanopy(klon)
       real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
       real petBcoef(klon), peqBcoef(klon)
-      real precip_rain(klon), precip_snow(klon)
       real lwdown(klon), swnet(klon), swdown(klon), ps(klon)
       real p1lay(klon)
@@ -545,10 +542,29 @@
 
 ! Parametres de sortie
-      real evap(klon), fluxsens(klon), fluxlat(klon)
+      real fluxsens(klon), fluxlat(klon)
       real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
       real emis_new(klon), z0_new(klon)
-      real dflux_l(klon), dflux_s(klon)
       real pctsrf_new(klon,nbsrf)
 c
+
+      if (.not. contreg) then
+        do k = 2, klev
+          do i = 1, knon
+            gamq(i,k) = 0.0
+            gamt(i,k) = 0.0
+          enddo
+        enddo
+      else
+        do k = 3, klev
+          do i = 1, knon
+            gamq(i,k)= 0.0
+            gamt(i,k)=  -1.0e-03
+          enddo
+        enddo
+        do i = 1, knon
+          gamq(i,2) = 0.0
+          gamt(i,2) = -2.5e-03
+        enddo
+      endif
 
       DO i = 1, knon
@@ -650,19 +666,28 @@
       spechum=q(:,1)
       p1lay = pplay(:,1)
+      zlev1 = delp(:,1)
+      swnet = swdown * (1 - albedo)
+c En attendant mieux
+      hum_air = 0.
+      ccanopy = 0.
+      tq_cdrag = 0.
 
       CALL interfsurf(itime, dtime, jour,
-     . klon, nisurf, knon, knindex, rlon, rlat,
+     . klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat,
      . debut, lafin, ok_veget,
-     . zlev, zlflu, u1lay, v1lay, temp_air, spechum, hum_air, ccanopy, 
+     . zlev1,  u1lay, v1lay, temp_air, spechum, hum_air, ccanopy, 
      . tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
      . precip_rain, precip_snow, lwdown, swnet, swdown,
-     . ts, p1lay, cal, beta, coef1lay, psref, radsol, dif_grnd,
+     . albedo, snow, qsol,
+     . ts, p1lay, coef1lay, psref, radsol,
      . ocean,
      . evap, fluxsens, fluxlat, dflux_l, dflux_s,              
-     . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
+     . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new, 
+     . zmasq)
 
       flux_t(:,1) = fluxsens
-      flux_q(:,1) = fluxlat
+      flux_q(:,1) = - evap
       d_ts = tsurf_new - ts
+      albedo = alb_new
 
 c==== une fois on a zx_h_ts, on peut faire l'iteration ========
@@ -1228,5 +1253,5 @@
       RETURN
       END
-      SUBROUTINE calbeta(dtime,indice,snow,qsol,
+      SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
      .                    vbeta,vcal,vdif)
       IMPLICIT none
@@ -1238,6 +1263,6 @@
 c Calculer quelques parametres pour appliquer la couche limite
 c ------------------------------------------------------------
-#include "dimensions.h"
-#include "dimphy.h"
+!#include "dimensions.h"
+!#include "dimphy.h"
 #include "YOMCST.h"
 #include "indicesol.h"
@@ -1256,13 +1281,13 @@
 c
       REAL dtime
-      REAL snow(klon,nbsrf), qsol(klon,nbsrf)
-      INTEGER indice
-C
-      REAL vbeta(klon)
-      REAL vcal(klon)
-      REAL vdif(klon)
+      REAL snow(knon,nbsrf), qsol(knon,nbsrf)
+      INTEGER indice, knon
+C
+      REAL vbeta(knon)
+      REAL vcal(knon)
+      REAL vdif(knon)
 C
       IF (indice.EQ.is_oce) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i)   = 0.0
           vbeta(i)  = 1.0
@@ -1272,5 +1297,5 @@
 c
       IF (indice.EQ.is_sic) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i) = calice
           IF (snow(i,is_sic) .GT. 0.0) vcal(i) = calsno
@@ -1282,5 +1307,5 @@
 c
       IF (indice.EQ.is_ter) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i) = calsol
           IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
@@ -1291,5 +1316,5 @@
 c
       IF (indice.EQ.is_lic) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i) = calice
           IF (snow(i,is_lic) .GT. 0.0) vcal(i) = calsno
