Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/pbl_surface_mod.F90	(revision 3807)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/pbl_surface_mod.F90	(revision 3823)
@@ -23,4 +23,7 @@
   USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
   USE coef_diff_turb_mod,  ONLY : coef_diff_turb
+  USE ioipsl_getin_p_mod,  ONLY : getin_p
+  USE cdrag_mod
+  USE stdlevvar_mod
   USE wx_pbl_mod,          ONLY : wx_pbl_init, wx_pbl_final, &
 !!                                  wx_pbl_fuse_no_dts, wx_pbl_split_no_dts, &
@@ -42,4 +45,6 @@
   INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
   !$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
+  integer, save :: iflag_new_t2mq2m
+  !$OMP THREADPRIVATE(iflag_new_t2mq2m)
 !FC
 !  integer, save :: iflag_frein
@@ -177,5 +182,5 @@
        alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
        alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
-       zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
+       zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout, &
        d_t,       d_q,       d_u,      d_v, d_t_diss, & 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -394,4 +399,5 @@
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
+    INTEGER, DIMENSION(klon, 6),  INTENT(OUT)       :: zn2mout ! nb of times t2m is out of the [tsol,temp]
     REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature 
@@ -537,4 +543,6 @@
     REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
     REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
+    INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
+    INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout, n2mout_x, n2mout_w
     REAL, DIMENSION(klon)              :: yustar
     REAL, DIMENSION(klon)              :: ywstar
@@ -838,4 +846,9 @@
 
     IF (first_call) THEN
+
+       iflag_new_t2mq2m=0
+       CALL getin_p('iflag_new_t2mq2m',iflag_new_t2mq2m)
+       print*,'pbl_iflag_new_t2mq2m=',iflag_new_t2mq2m
+
        print*,'PBL SURFACE AVEC GUSTINESS'
        first_call=.FALSE.
@@ -903,4 +916,5 @@
  zxfluxlat(:)=0.
  zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
+ zn2mout(:,:)=0 ;
  d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
  zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
@@ -1821,9 +1835,15 @@
 
           ! Calculate the temperature et relative humidity at 2m and the wind at 10m 
+          IF (iflag_new_t2mq2m==1) THEN
+           CALL stdlevvarn(klon, knon, is_ter, zxli, &
+               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
+               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
+               yt2m, yq2m, yt10m, yq10m, yu10m, yn2mout(:, nsrf, :))
+          ELSE
           CALL stdlevvar(klon, knon, is_ter, zxli, &
                yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
                yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
                yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
-          
+          ENDIF
        ENDIF
 
@@ -2614,9 +2634,26 @@
 !!! jyg le 07/02/2012
        IF (iflag_split .eq.0) THEN
+        IF (iflag_new_t2mq2m==1) THEN
+         CALL stdlevvarn(klon, knon, nsrf, zxli, &
+            uzon, vmer, tair1, qair1, zgeo1, &
+            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
+            yt2m, yq2m, yt10m, yq10m, yu10m, yn2mout(:, nsrf, :))
+        ELSE
         CALL stdlevvar(klon, knon, nsrf, zxli, &
             uzon, vmer, tair1, qair1, zgeo1, &
             tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
             yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
+        ENDIF
        ELSE  !(iflag_split .eq.0)
+        IF (iflag_new_t2mq2m==1) THEN
+         CALL stdlevvarn(klon, knon, nsrf, zxli, &
+            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
+            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
+            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yn2mout_x(:, nsrf, :))
+         CALL stdlevvarn(klon, knon, nsrf, zxli, &
+            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
+            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
+            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yn2mout_w(:, nsrf, :))
+        ELSE
         CALL stdlevvar(klon, knon, nsrf, zxli, &
             uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
@@ -2627,4 +2664,5 @@
             tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
             yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
+        ENDIF
 !!!
        ENDIF  ! (iflag_split .eq.0)
@@ -2640,4 +2678,9 @@
           u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
           v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
+!
+          DO k = 1, 6
+           n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
+          END DO
+!
         ENDDO
        ELSE  !(iflag_split .eq.0)
@@ -2650,4 +2693,9 @@
           u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
           v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
+!
+          DO k = 1, 6
+           n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
+          END DO
+!
         ENDDO
         DO j=1, knon
@@ -2663,4 +2711,9 @@
           u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
           v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
+!
+          DO k = 1, 6
+           n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
+          END DO
+!
         ENDDO
 !!!
@@ -2919,5 +2972,5 @@
 !
     zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
-    zt2m(:) = 0.0    ; zq2m(:) = 0.0 
+    zt2m(:) = 0.0    ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
     zustar(:)=0.0 ; zu10m(:) = 0.0   ; zv10m(:) = 0.0
     s_pblh(:) = 0.0  ; s_plcl(:) = 0.0 
@@ -2972,4 +3025,9 @@
           zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
           zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
+!
+          DO k = 1, 6
+           zn2mout(i,k)  = zn2mout(i,k)  + n2mout(i,nsrf,k)  * pctsrf(i,nsrf)
+          ENDDO
+!
           zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
           wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
@@ -3003,4 +3061,9 @@
           zt2m(i)  = zt2m(i)  + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)
           zq2m(i)  = zq2m(i)  + q2m_x(i,nsrf)  * pctsrf(i,nsrf)
+!
+          DO k = 1, 6
+           zn2mout(i,k)  = zn2mout(i,k)  + n2mout_x(i,nsrf,k)  * pctsrf(i,nsrf)
+          ENDDO
+!
           zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
           wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
@@ -3081,5 +3144,5 @@
     DO nsrf = 1, nbsrf
        DO i = 1, klon
-          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
+          zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
           zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
        ENDDO
