Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5015)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5022)
@@ -33,4 +33,7 @@
                                   wx_pbl_check, wx_pbl_dts_check, wx_evappot
   use config_ocean_skin_m, only: activate_ocean_skin
+#ifdef ISO
+  USE infotrac_phy, ONLY: niso,ntraciso=>ntiso    
+#endif
 
   IMPLICIT NONE
@@ -49,4 +52,13 @@
   !$OMP THREADPRIVATE(ydTs0, ydqs0)
 
+#ifdef ISO
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: xtsnow   ! snow at surface
+  !$OMP THREADPRIVATE(xtsnow)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: Rland_ice   ! snow at surface
+  !$OMP THREADPRIVATE(Rland_ice)  
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: Roce   ! snow at surface
+  !$OMP THREADPRIVATE(Roce)  
+#endif
+
   INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
   !$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
@@ -178,4 +190,69 @@
 
   END SUBROUTINE pbl_surface_init
+
+#ifdef ISO
+  SUBROUTINE pbl_surface_init_iso(xtsnow_rst,Rland_ice_rst)
+
+! This routine should be called after the restart file has been read.
+! This routine initialize the restart variables and does some validation tests
+! for the index of the different surfaces and tests the choice of type of ocean.
+
+    USE indice_sol_mod
+    USE print_control_mod, ONLY: lunout
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    USE isotopes_verif_mod
+#endif
+    IMPLICIT NONE
+
+    INCLUDE "dimsoil.h"
+ 
+! Input variables
+!****************************************************************************************
+    REAL, DIMENSION(niso,klon, nbsrf), INTENT(IN)          :: xtsnow_rst
+    REAL, DIMENSION(niso,klon), INTENT(IN)          :: Rland_ice_rst
+  
+! Local variables
+!****************************************************************************************
+    INTEGER                       :: ierr
+    CHARACTER(len=80)             :: abort_message
+    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
+    integer i,ixt
+    
+!****************************************************************************************
+! Allocate and initialize module variables with fields read from restart file.
+!
+!****************************************************************************************    
+
+    ALLOCATE(xtsnow(niso,klon,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(Rland_ice(niso,klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(Roce(niso,klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    xtsnow(:,:,:)  = xtsnow_rst(:,:,:)
+    Rland_ice(:,:) = Rland_ice_rst(:,:)
+    Roce(:,:)      = 0.0
+
+#ifdef ISOVERIF 
+      IF (iso_eau >= 0) THEN
+         CALL iso_verif_egalite_vect2D( &
+     &           xtsnow,snow, &
+     &           'pbl_surface_mod 170',niso,klon,nbsrf)
+         DO i=1,klon  
+            IF (iso_eau >= 0) THEN  
+              CALL iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
+     &         'pbl_surf_mod 177')
+            ENDIF
+         ENDDO
+      ENDIF
+#endif
+
+  END SUBROUTINE pbl_surface_init_iso
+#endif
+
 !  
 !****************************************************************************************
@@ -241,5 +318,12 @@
 !FC
 !!!
-                        )
+#ifdef ISO
+     &   ,xtrain_f, xtsnow_f,xt, &
+     &   wake_dlxt,zxxtevap,xtevap, &
+     &   d_xt,d_xt_w,d_xt_x, &
+     &   xtsol,dflux_xt,zxxtsnow,zxfluxxt,flux_xt, &
+     &   h1_diag,runoff_diag,xtrunoff_diag &
+#endif      
+     &   )
 !****************************************************************************************
 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
@@ -316,4 +400,13 @@
     USE mod_grid_phy_lmdz,  ONLY : nbp_lon, nbp_lat, grid1dto2d_glo
     USE print_control_mod,  ONLY : prt_level,lunout
+#ifdef ISO
+  USE isotopes_mod, ONLY: Rdefault,iso_eau
+#ifdef ISOVERIF
+        USE isotopes_verif_mod
+#endif
+#ifdef ISOTRAC
+        USE isotrac_mod, only: index_iso
+#endif
+#endif
     USE ioipsl_getin_p_mod, ONLY : getin_p
     use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &
@@ -368,5 +461,11 @@
     REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
 
-    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud fraction
+    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud 
+
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: xt       ! water vapour (kg/kg)
+    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtrain_f  ! rain fall
+    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtsnow_f  ! snow fall
+#endif
 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -381,5 +480,7 @@
     REAL, DIMENSION(klon),        INTENT(IN)        :: wake_dens
 !!!
-
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: wake_dlxt   
+#endif
 ! Input/Output variables
 !****************************************************************************************
@@ -450,4 +551,14 @@
     REAL, INTENT(OUT):: zcoefm(:, :, :) ! (klon, klev, nbsrf + 1)
     ! coef for turbulent diffusion of U and V (?), mean for each grid point
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon),        INTENT(OUT)       :: zxxtevap     ! water vapour flux at surface, positiv upwards
+    REAL, DIMENSION(ntraciso,klon, klev),  INTENT(OUT)       :: d_xt        ! change in water vapour
+    REAL, DIMENSION(klon),                 INTENT(OUT)       :: runoff_diag
+    REAL, DIMENSION(niso,klon),            INTENT(OUT)       :: xtrunoff_diag
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_w
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_x
+#endif
+
+
 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -513,5 +624,10 @@
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
 !FC
-    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT)  :: treedrg      ! tree drag (m)               
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg  ! tree drag (m)               
+#ifdef ISO        
+    REAL, DIMENSION(niso,klon),   INTENT(OUT)       :: xtsol      ! water height in the soil (mm)
+    REAL, DIMENSION(ntraciso,klon, nbsrf)           :: xtevap     ! evaporation at surface
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: h1_diag    ! just diagnostic, not useful
+#endif
 
 
@@ -527,4 +643,10 @@
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs   ! blowind snow vertical flux (kg/m**2
 
+#ifdef ISO   
+    REAL, DIMENSION(ntraciso,klon),              INTENT(OUT) :: dflux_xt    ! change of water vapour flux
+    REAL, DIMENSION(niso,klon),                  INTENT(OUT) :: zxxtsnow    ! snow at surface, mean for each grid point
+    REAL, DIMENSION(ntraciso,klon, klev),        INTENT(OUT) :: zxfluxxt    ! water vapour flux, mean for each grid point 
+    REAL, DIMENSION(ntraciso,klon, klev, nbsrf), INTENT(OUT) :: flux_xt     ! water vapour flux(latent flux) (kg/m**2/s)  
+#endif
 
 ! Martin
@@ -575,4 +697,12 @@
     REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
     REAL, DIMENSION(klon)              :: yrain_f, ysnow_f, ybs_f
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)     :: yxt1
+    REAL, DIMENSION(niso,klon)         :: yxtsnow, yxtsol   
+    REAL, DIMENSION(ntraciso,klon)     :: yxtrain_f, yxtsnow_f 
+    REAL, DIMENSION(klon)              :: yrunoff_diag
+    REAL, DIMENSION(niso,klon)         :: yxtrunoff_diag
+    REAL, DIMENSION(niso,klon)         :: yRland_ice    
+#endif
     REAL, DIMENSION(klon)              :: ysolsw, ysollw
     REAL, DIMENSION(klon)              :: yfder
@@ -583,4 +713,8 @@
     REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
     REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)     ::  y_flux_xt1
+    REAL, DIMENSION(ntraciso,klon)     ::  y_dflux_xt
+#endif
     REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
     REAL, DIMENSION(klon)              :: y_flux_bs, y_flux0
@@ -610,4 +744,7 @@
     REAL, DIMENSION(klon)              :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
     REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)     :: AcoefXT, BcoefXT
+#endif
     REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
     REAL, DIMENSION(klon)              :: AcoefQBS, BcoefQBS
@@ -628,4 +765,10 @@
     REAL, DIMENSION(klon,klev)         :: yu, yv
     REAL, DIMENSION(klon,klev)         :: yt, yq, yqbs
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)      :: yxtevap
+    REAL, DIMENSION(ntraciso,klon,klev) :: y_d_xt
+    REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt
+    REAL, DIMENSION(ntraciso,klon,klev) :: yxt   
+#endif
     REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
     REAL, DIMENSION(klon,klev)         :: delp
@@ -699,4 +842,17 @@
     REAL, DIMENSION(klon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
     REAL, DIMENSION(klon)              :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,klev)         :: yxt_x, yxt_w
+    REAL, DIMENSION(ntraciso,klon)              :: y_flux_xt1_x , y_flux_xt1_w   
+    REAL, DIMENSION(ntraciso,klon,klev)         :: y_flux_xt_x,y_d_xt_x,zxfluxxt_x
+    REAL, DIMENSION(ntraciso,klon,klev)         :: y_flux_xt_w,y_d_xt_w,zxfluxxt_w
+    REAL, DIMENSION(ntraciso,klon,klev,nbsrf)   :: flux_xt_x, flux_xt_w
+    REAL, DIMENSION(ntraciso,klon)              :: AcoefXT_x, BcoefXT_x
+    REAL, DIMENSION(ntraciso,klon)              :: AcoefXT_w, BcoefXT_w
+    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT, DcoefXT
+    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT_x, DcoefXT_x
+    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT_w, DcoefXT_w
+    REAL, DIMENSION(ntraciso,klon,klev)         :: gama_xt,gama_xt_x,gama_xt_w
+#endif
 !!!
 !!!jyg le 08/02/2012
@@ -891,5 +1047,5 @@
     REAL, DIMENSION(klon)              :: yrmu0
     ! Martin
-    REAL, DIMENSIOn(klon)              :: yri0
+    REAL, DIMENSION(klon)              :: yri0
 
     REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
@@ -898,4 +1054,12 @@
     ! dt_ds, tkt, tks, taur, sss on ocean points
     REAL :: missing_val
+#ifdef ISO
+    REAL, DIMENSION(klon)       :: h1
+    INTEGER                     :: ixt
+!#ifdef ISOVERIF
+!    integer iso_verif_positif_nostop
+!#endif    
+#endif
+
 !****************************************************************************************
 ! End of declarations
@@ -925,4 +1089,42 @@
       iflag_split_ref = mod(iflag_pbl_split,10)
       iflag_split = iflag_split_ref
+
+#ifdef ISO      
+#ifdef ISOVERIF
+      DO i=1,klon
+        DO ixt=1,niso
+          CALL iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 608')
+        ENDDO
+      ENDDO
+#endif
+#ifdef ISOVERIF
+      DO i=1,klon  
+        IF (iso_eau >= 0) THEN  
+          CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
+     &         'pbl_surf_mod 585',errmax,errmaxrel)
+          CALL iso_verif_egalite_choix(xtsnow_f(iso_eau,i),snow_f(i), &
+     &         'pbl_surf_mod 594',errmax,errmaxrel)
+          IF (iso_verif_egalite_choix_nostop(xtsol(iso_eau,i),qsol(i), &
+     &         'pbl_surf_mod 596',errmax,errmaxrel) == 1) THEN
+                WRITE(*,*) 'i=',i
+                STOP
+          ENDIF
+          DO nsrf=1,nbsrf
+            CALL iso_verif_egalite_choix(xtsnow(iso_eau,i,nsrf),snow(i,nsrf), &
+     &         'pbl_surf_mod 598',errmax,errmaxrel)
+          ENDDO
+        ENDIF !IF (iso_eau >= 0) THEN   
+      ENDDO !DO i=1,knon  
+      DO k=1,klev
+        DO i=1,klon  
+          IF (iso_eau >= 0) THEN  
+            CALL iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &
+     &           'pbl_surf_mod 595',errmax,errmaxrel)
+          ENDIF !IF (iso_eau >= 0) THEN  
+        ENDDO !DO i=1,knon  
+      ENDDO !DO k=1,klev
+#endif
+#endif
+
 
 !****************************************************************************************
@@ -990,4 +1192,14 @@
       PRINT*,'WARNING : On impose qsol=',qsol0
       qsol(:)=qsol0
+#ifdef ISO
+      DO ixt=1,niso
+        xtsol(ixt,:)=qsol0*Rdefault(ixt)
+      ENDDO
+#ifdef ISOTRAC      
+      DO ixt=1+niso,ntraciso
+        xtsol(ixt,:)=qsol0*Rdefault(index_iso(ixt))
+      ENDDO
+#endif       
+#endif
     ENDIF
 !****************************************************************************************
@@ -1040,4 +1252,13 @@
  qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
  runoff(:)=0.
+#ifdef ISO
+zxxtevap(:,:)=0.
+ d_xt(:,:,:)=0. 
+ d_xt_x(:,:,:)=0.
+ d_xt_w(:,:,:)=0.
+ flux_xt(:,:,:,:)=0. 
+! xtsnow(:,:,:)=0.! attention, xtsnow est l'équivalent de snow et non de qsnow
+ xtevap(:,:,:)=0.
+#endif
     IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
        zcoefh(:,:,:) = 0.0
@@ -1129,4 +1350,17 @@
 !FC
 
+#ifdef ISO
+   yxtrain_f = 0.0 ; yxtsnow_f = 0.0
+   yxtsnow  = 0.0
+   yxt = 0.0
+   yxtsol = 0.0
+   flux_xt = 0.0
+   yRland_ice = 0.0
+!   d_xt = 0.0      
+   y_dflux_xt = 0.0  
+   dflux_xt=0.0 
+   y_d_xt_x=0.      ; y_d_xt_w=0.       
+#endif 
+
 ! >> PC
 !the yfields_out variable is defined in (klon,nbcf_out) even if it is used on
@@ -1155,4 +1389,7 @@
     fluxlat_x(:,:)=0. ;           fluxlat_w(:,:)=0. 
 !>jyg
+#ifdef ISO
+    flux_xt_x(:,:,:,:)=0. ;          flux_xt_w(:,:,:,:)=0.
+#endif
 !
 !jyg<
@@ -1454,4 +1691,34 @@
           yfluxbs(j)=0.0
           y_flux_bs(j) = 0.0
+!!!
+#ifdef ISO
+          DO ixt=1,ntraciso
+            yxtrain_f(ixt,j) = xtrain_f(ixt,i)
+            yxtsnow_f(ixt,j) = xtsnow_f(ixt,i)  
+          ENDDO
+          DO ixt=1,niso
+            yxtsnow(ixt,j)   = xtsnow(ixt,i,nsrf)
+          ENDDO    
+          !IF (nsrf == is_lic) THEN
+          DO ixt=1,niso
+            yRland_ice(ixt,j)= Rland_ice(ixt,i)  
+          ENDDO    
+          !endif !IF (nsrf == is_lic) THEN
+#ifdef ISOVERIF
+          IF (iso_eau >= 0) THEN
+              call iso_verif_egalite_choix(ysnow_f(j), &
+     &          yxtsnow_f(iso_eau,j),'pbl_surf_mod 862', &
+     &          errmax,errmaxrel)
+              call iso_verif_egalite_choix(ysnow(j), &
+     &          yxtsnow(iso_eau,j),'pbl_surf_mod 872', &
+     &          errmax,errmaxrel)
+          ENDIF
+#endif
+#ifdef ISOVERIF
+         DO ixt=1,ntraciso
+           call iso_verif_noNaN(yxtsnow_f(ixt,j),'pbl_surf_mod 921')
+         ENDDO
+#endif
+#endif
        ENDDO
 ! >> PC
@@ -1493,4 +1760,9 @@
              yq(j,k) = q(i,k)
              yqbs(j,k)=qbs(i,k)
+#ifdef ISO
+             DO ixt=1,ntraciso   
+               yxt(ixt,j,k) = xt(ixt,i,k)
+             ENDDO !DO ixt=1,ntraciso
+#endif
           ENDDO
         ENDDO
@@ -1510,4 +1782,10 @@
              yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
 !!!
+#ifdef ISO
+             DO ixt=1,ntraciso
+               yxt_x(ixt,j,k) = xt(ixt,i,k)-wake_s(i)*wake_dlxt(ixt,i,k)
+               yxt_w(ixt,j,k) = xt(ixt,i,k)+(1.-wake_s(i))*wake_dlxt(ixt,i,k)
+             ENDDO
+#endif
           ENDDO
         ENDDO
@@ -1565,4 +1843,9 @@
              i = ni(j)
              yqsol(j) = qsol(i)
+#ifdef ISO
+             DO ixt=1,niso
+               yxtsol(ixt,j) = xtsol(ixt,i)
+             ENDDO
+#endif
           ENDDO
        ENDIF
@@ -1859,5 +2142,9 @@
             Kcoef_hq, gama_q, gama_h, &
 !!!
-            AcoefH, AcoefQ, BcoefH, BcoefQ)
+            AcoefH, AcoefQ, BcoefH, BcoefQ &
+#ifdef ISO
+         &   ,yxt, CcoefXT, DcoefXT, gama_xt, AcoefXT, BcoefXT & 
+#endif               
+         &   )
        ELSE  !(iflag_split .eq.0)
         CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
@@ -1867,5 +2154,9 @@
             Kcoef_hq_x, gama_q_x, gama_h_x, &
 !!!
-            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
+            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x &
+#ifdef ISO
+         &   ,yxt_x, CcoefXT_x, DcoefXT_x, gama_xt_x, AcoefXT_x, BcoefXT_x & 
+#endif               
+         &   )
 !!!
        IF (prt_level >=10) THEN
@@ -1882,5 +2173,9 @@
             Kcoef_hq_w, gama_q_w, gama_h_w, &
 !!!
-            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
+            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w &
+#ifdef ISO
+         &   ,yxt_w, CcoefXT_w, DcoefXT_w, gama_xt_w, AcoefXT_w, BcoefXT_w & 
+#endif               
+         &   )
 !!!
        IF (prt_level >=10) THEN
@@ -1964,5 +2259,13 @@
          yt1(:) = yt(:,1)
          yq1(:) = yq(:,1)
+#ifdef ISO
+         yxt1(:,:) = yxt(:,:,1)
+#endif
+
        ELSE IF (iflag_split .ge. 1) THEN
+#ifdef ISO
+        call abort_gcm('pbl_surface_mod 2149','isos pas encore dans iflag_split=1',1)
+#endif
+
 !
 ! Cdragq computation
@@ -2126,5 +2429,11 @@
                yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
                y_flux_u1, y_flux_v1, &
-               yveget,ylai,yheight )
+               yveget,ylai,yheight   &
+#ifdef ISO
+         &      ,yxtrain_f, yxtsnow_f,yxt1, &
+         &      yxtsnow,yxtsol,yxtevap,h1, &
+         &      yrunoff_diag,yxtrunoff_diag,yRland_ice &
+#endif               
+         &      )
  
 !FC quid qd yveget ylai yheight ne sont pas definit
@@ -2156,5 +2465,27 @@
           ENDDO
       ENDIF
-     
+
+#ifdef ISOVERIF
+        DO j=1,knon
+          DO ixt=1,ntraciso
+            CALL iso_verif_noNaN(yxtevap(ixt,j), &
+         &      'pbl_surface 1056a: apres surf_land')
+          ENDDO
+          DO ixt=1,niso
+            CALL iso_verif_noNaN(yxtsol(ixt,j), &
+         &      'pbl_surface 1056b: apres surf_land')
+          ENDDO
+        ENDDO
+#endif
+#ifdef ISOVERIF
+!        write(*,*) 'pbl_surface_mod 1038: sortie surf_land'
+        DO j=1,knon
+          IF (iso_eau >= 0) THEN     
+                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
+     &                                  ysnow(j),'pbl_surf_mod 1043')
+          ENDIF !if (iso_eau.gt.0) then
+        ENDDO !DO i=1,klon
+#endif
+    
        CASE(is_lic)
           ! Martin
@@ -2177,5 +2508,10 @@
                   ysnowhgt, yqsnow, ytoice, ysissnow, &
                   yalb3_new, yrunoff, &
-                  y_flux_u1, y_flux_v1)
+                  y_flux_u1, y_flux_v1 &
+#ifdef ISO
+                  &    ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice &
+                  &    ,yxtsnow,yxtsol,yxtevap &
+#endif              
+                  &    )
              
              !jyg<
@@ -2199,5 +2535,27 @@
                 ENDDO
              ENDIF
-             
+
+#ifdef ISOVERIF
+             DO j=1,knon
+               DO ixt=1,ntraciso
+                 CALL iso_verif_noNaN(yxtevap(ixt,j), &
+                        &             'pbl_surface 1095a: apres surf_landice')
+               ENDDO
+                do ixt=1,niso
+                   call iso_verif_noNaN(yxtsol(ixt,j), &
+                        &      'pbl_surface 1095b: apres surf_landice')
+                enddo
+             enddo
+#endif
+#ifdef ISOVERIF
+             !write(*,*) 'pbl_surface_mod 1060: sortie surf_landice'
+             do j=1,knon
+               IF (iso_eau >= 0) THEN     
+                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
+                        &               ysnow(j),'pbl_surf_mod 1064')
+               ENDIF !if (iso_eau >= 0) THEN
+             ENDDO !DO i=1,klon
+#endif
+            
           END IF
           
@@ -2216,5 +2574,10 @@
                y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
                yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
-               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
+               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss &
+#ifdef ISO
+         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
+         &      yxtsnow,yxtevap,h1 &
+#endif               
+         &      )
       IF (prt_level >=10) THEN
           print *,'arg de surf_ocean: ycdragh ',ycdragh(1:knon)
@@ -2257,5 +2620,10 @@
 !albedo SB <<<
                ytsurf_new, y_dflux_t, y_dflux_q, &
-               y_flux_u1, y_flux_v1)
+               y_flux_u1, y_flux_v1 &
+#ifdef ISO
+         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
+         &      yxtsnow,yxtsol,yxtevap,Rland_ice &
+#endif               
+         &      )
           
 ! Special DICE MPL 05082013 puis BOMEX MPL 20150410
@@ -2265,5 +2633,27 @@
           y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
           ENDDO
-      ENDIF
+       ENDIF
+
+#ifdef ISOVERIF
+        DO j=1,knon
+          DO ixt=1,ntraciso
+            CALL iso_verif_noNaN(yxtevap(ixt,j), &
+         &                       'pbl_surface 1165a: apres surf_seaice')
+          ENDDO
+          DO ixt=1,niso
+            CALL iso_verif_noNaN(yxtsol(ixt,j), &
+         &      'pbl_surface 1165b: apres surf_seaice')
+          ENDDO
+        ENDDO
+#endif
+#ifdef ISOVERIF
+        !write(*,*) 'pbl_surface_mod 1077: sortie surf_seaice'
+        DO j=1,knon
+          IF (iso_eau >= 0) THEN     
+                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
+     &                                  ysnow(j),'pbl_surf_mod 1106')
+          ENDIF !IF (iso_eau >= 0) THEN
+        ENDDO !DO i=1,klon
+#endif
 
        CASE DEFAULT
@@ -2335,4 +2725,7 @@
           y_flux_t1(j) =  yfluxsens(j)
           y_flux_q1(j) = -yevap(j)
+#ifdef ISO
+          y_flux_xt1(:,:) = -yxtevap(:,:)
+#endif
           ENDDO
         ENDIF ! (ok_flux_surf)
@@ -2350,4 +2743,8 @@
 
        IF (iflag_split .GE. 1) THEN
+#ifdef ISO
+        call abort_gcm('pbl_surface_mod 2607','isos pas encore dans iflag_split=1',1)
+#endif
+!
 !
          IF (nsrf .ne. is_oce) THEN
@@ -2567,5 +2964,11 @@
             Kcoef_hq, gama_q, gama_h, &
 !!!
-            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))    
+            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:) &
+#ifdef ISO
+        &    ,yxt,y_flux_xt1 &
+        &    ,AcoefXT,BcoefXT,CcoefXT,DcoefXT,gama_xt &
+        &    ,y_flux_xt(:,:,:),y_d_xt(:,:,:) &
+#endif
+        &    )    
        ELSE  !(iflag_split .eq.0)
         CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
@@ -2576,5 +2979,11 @@
             Kcoef_hq_x, gama_q_x, gama_h_x, &
 !!!
-            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))    
+            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:) &
+#ifdef ISO
+        &    ,yxt_x,y_flux_xt1_x &
+        &    ,AcoefXT_x,BcoefXT_x,CcoefXT_x,DcoefXT_x,gama_xt_x &
+        &    ,y_flux_xt_x(:,:,:),y_d_xt_x(:,:,:) &
+#endif
+        &    )    
 !
        CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
@@ -2585,5 +2994,11 @@
             Kcoef_hq_w, gama_q_w, gama_h_w, &
 !!!
-            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))    
+            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:) &
+#ifdef ISO
+        &    ,yxt_w,y_flux_xt1_w &
+        &    ,AcoefXT_w,BcoefXT_w,CcoefXT_w,DcoefXT_w,gama_xt_w &
+        &    ,y_flux_xt_w(:,:,:),y_d_xt_w(:,:,:) &
+#endif
+        &    )    
 !!!
        ENDIF  ! (iflag_split .eq.0)
@@ -2703,6 +3118,25 @@
              flux_u(i,k,nsrf) = y_flux_u(j,k)
              flux_v(i,k,nsrf) = y_flux_v(j,k)
+
+#ifdef ISO
+             DO ixt=1,ntraciso
+                y_d_xt(ixt,j,k)  = y_d_xt(ixt,j,k) * ypct(j)
+                flux_xt(ixt,i,k,nsrf) = y_flux_xt(ixt,j,k)
+             ENDDO ! DO ixt=1,ntraciso
+             h1_diag(i)=h1(j)
+#endif
+
            ENDDO
         ENDDO
+
+#ifdef ISO
+#ifdef ISOVERIF
+        if (iso_eau.gt.0) then
+         call iso_verif_egalite_vect2D( &
+                y_d_xt,y_d_q, &
+                'pbl_surface_mod 2600',ntraciso,klon,klev)
+        endif       
+#endif
+#endif
 
        ELSE  !(iflag_split .eq.0)
@@ -2722,4 +3156,11 @@
             flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
             flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
+
+#ifdef ISO
+            DO ixt=1,ntraciso
+              y_d_xt_x(ixt,j,k)  = y_d_xt_x(ixt,j,k) * ypct(j)
+              flux_xt_x(ixt,i,k,nsrf) = y_flux_xt_x(ixt,j,k)
+            ENDDO ! DO ixt=1,ntraciso
+#endif
           ENDDO
         ENDDO
@@ -2739,4 +3180,12 @@
             flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
             flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
+
+#ifdef ISO
+            DO ixt=1,ntraciso
+              y_d_xt_w(ixt,j,k)  = y_d_xt_w(ixt,j,k) * ypct(j)
+              flux_xt_w(ixt,i,k,nsrf) = y_flux_xt_w(ixt,j,k)
+            ENDDO ! do ixt=1,ntraciso
+#endif
+
           ENDDO
         ENDDO
@@ -2750,4 +3199,9 @@
             flux_u(i,k,nsrf) = flux_u_x(i,k,nsrf)+ywake_s(j)*(flux_u_w(i,k,nsrf)-flux_u_x(i,k,nsrf))
             flux_v(i,k,nsrf) = flux_v_x(i,k,nsrf)+ywake_s(j)*(flux_v_w(i,k,nsrf)-flux_v_x(i,k,nsrf))
+#ifdef ISO
+            DO ixt=1,ntraciso
+              flux_xt(ixt,i,k,nsrf) = flux_xt_x(ixt,i,k,nsrf)+ywake_s(j)*(flux_xt_w(ixt,i,k,nsrf)-flux_xt_x(ixt,i,k,nsrf))
+            ENDDO ! do ixt=1,ntraciso
+#endif
           ENDDO
         ENDDO
@@ -2807,4 +3261,24 @@
           dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
           dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
+#ifdef ISO
+        DO ixt=1,niso
+          xtsnow(ixt,i,nsrf) = yxtsnow(ixt,j)  
+        ENDDO
+        DO ixt=1,ntraciso
+          xtevap(ixt,i,nsrf) = - flux_xt(ixt,i,1,nsrf)
+          dflux_xt(ixt,i) = dflux_xt(ixt,i) + y_dflux_xt(ixt,j)*ypct(j)
+        ENDDO  
+        IF (nsrf == is_lic) THEN
+          DO ixt=1,niso
+            Rland_ice(ixt,i) = yRland_ice(ixt,j)  
+          ENDDO
+        ENDIF !IF (nsrf == is_lic) THEN     
+#ifdef ISOVERIF
+        IF (iso_eau.gt.0) THEN  
+          call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
+     &         'pbl_surf_mod 1230',errmax,errmaxrel)
+        ENDIF !if (iso_eau.gt.0) then
+#endif        
+#endif
        ENDDO
 
@@ -2911,4 +3385,11 @@
              i = ni(j)
              qsol(i) = yqsol(j)
+#ifdef ISO
+             runoff_diag(i)=yrunoff_diag(j)   
+             DO ixt=1,niso
+               xtsol(ixt,i) = yxtsol(ixt,j)
+               xtrunoff_diag(ixt,i)=yxtrunoff_diag(ixt,j)
+             ENDDO
+#endif
           ENDDO
        ENDIF
@@ -2923,5 +3404,22 @@
           ENDDO
        ENDDO
-       
+
+#ifdef ISO
+#ifdef ISOVERIF
+       !write(*,*) 'pbl_surface 2858'
+       DO i = 1, klon
+         DO ixt=1,niso
+           call iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 1405')
+         ENDDO
+       ENDDO
+#endif
+#ifdef ISOVERIF
+     IF (iso_eau.gt.0) THEN
+        call iso_verif_egalite_vect2D( &
+                y_d_xt,y_d_q, &
+                'pbl_surface_mod 1261',ntraciso,klon,klev)
+     ENDIF !if (iso_eau.gt.0) then
+#endif
+#endif
 !!! jyg le 07/02/2012
        IF (iflag_split .ge.1) THEN
@@ -2942,4 +3440,11 @@
            d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
            d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
+#ifdef ISO
+           DO ixt=1,ntraciso
+             d_xt_x(ixt,i,k) = d_xt_x(ixt,i,k) + y_d_xt_x(ixt,j,k)
+             d_xt_w(ixt,i,k) = d_xt_w(ixt,i,k) + y_d_xt_w(ixt,j,k)
+           ENDDO ! DO ixt=1,ntraciso
+#endif
+
 !
 !!           d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
@@ -2957,4 +3462,9 @@
              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)
+#ifdef ISO
+             DO ixt=1,ntraciso
+               d_xt(ixt,i,k) = d_xt(ixt,i,k) + y_d_xt(ixt,j,k)
+             ENDDO !DO ixt=1,ntraciso
+#endif
              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)
@@ -2971,4 +3481,24 @@
          ENDDO
         ENDIF
+
+#ifdef ISO
+#ifdef ISOVERIF
+!        write(*,*) 'd_q,d_xt(iso_eau,554,19)=',d_q(554,19),d_xt(iso_eau,554,19)
+!        write(*,*) 'pbl_surface 2929: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
+!        write(*,*) 'y_d_q,y_d_xt(iso_eau,2,1)=',y_d_q(2,1),y_d_xt(iso_eau,2,1)
+!        write(*,*) 'iso_eau.gt.0=',iso_eau.gt.0
+        call iso_verif_noNaN_vect2D( &
+     &           d_xt, &
+     &           'pbl_surface 1385',ntraciso,klon,klev)  
+     IF (iso_eau >= 0) THEN
+        call iso_verif_egalite_vect2D( &
+                y_d_xt,y_d_q, &
+                'pbl_surface_mod 2945',ntraciso,klon,klev)
+        call iso_verif_egalite_vect2D( &
+                d_xt,d_q, &
+                'pbl_surface_mod 1276',ntraciso,klon,klev)
+     ENDIF !IF (iso_eau >= 0) THEN
+#endif
+#endif
 
 !      print*,'Dans pbl OK4'
@@ -3358,4 +3888,15 @@
    iflag_split=iflag_split_ref
 
+#ifdef ISO
+#ifdef ISOVERIF
+!        write(*,*) 'pbl_surface tmp 3249: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
+    IF (iso_eau >= 0) THEN
+        call iso_verif_egalite_vect2D( &
+                d_xt,d_q, &
+                'pbl_surface_mod 1276',ntraciso,klon,klev)
+    ENDIF !IF (iso_eau >= 0) THEN
+#endif
+#endif
+
 !****************************************************************************************
 ! 16) Calculate the mean value over all sub-surfaces for some variables
@@ -3379,4 +3920,10 @@
     zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
     zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
+#ifdef ISO
+      zxfluxxt(:,:,:) = 0.0 
+      zxfluxxt_x(:,:,:) = 0.0
+      zxfluxxt_w(:,:,:) = 0.0
+#endif
+
 
 !!! jyg le 07/02/2012
@@ -3397,4 +3944,10 @@
               zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
               zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
+#ifdef ISO
+              DO ixt=1,ntraciso
+                zxfluxxt_x(ixt,i,k) = zxfluxxt_x(ixt,i,k) + flux_xt_x(ixt,i,k,nsrf) * pctsrf(i,nsrf)
+                zxfluxxt_w(ixt,i,k) = zxfluxxt_w(ixt,i,k) + flux_xt_w(ixt,i,k,nsrf) * pctsrf(i,nsrf)
+              ENDDO ! DO ixt=1,ntraciso
+#endif
             ENDDO
           ENDDO
@@ -3416,4 +3969,9 @@
              zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
+#ifdef ISO
+             DO ixt=1,niso
+               zxfluxxt(ixt,i,k) = zxfluxxt(ixt,i,k) + flux_xt(ixt,i,k,nsrf) * pctsrf(i,nsrf)
+             ENDDO ! DO ixt=1,niso
+#endif
           ENDDO
        ENDDO
@@ -3440,4 +3998,12 @@
        END DO
     endif
+
+#ifdef ISO
+    DO i = 1, klon
+      DO ixt=1,ntraciso
+        zxxtevap(ixt,i)     = - zxfluxxt(ixt,i,1)
+      ENDDO
+    ENDDO
+#endif
 
 !!!
@@ -3615,8 +4181,17 @@
     zxqsurf(:) = 0.0
     zxsnow(:)  = 0.0
+#ifdef ISO
+    zxxtsnow(:,:)  = 0.0
+#endif
+
     DO nsrf = 1, nbsrf
        DO i = 1, klon
           zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
           zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
+#ifdef ISO
+          DO ixt=1,niso
+            zxxtsnow(ixt,i)  = zxxtsnow(ixt,i)  + xtsnow(ixt,i,nsrf)  * pctsrf(i,nsrf)
+          ENDDO ! DO ixt=1,niso
+#endif
        ENDDO
     ENDDO
@@ -3630,7 +4205,17 @@
 !****************************************************************************************
 !
-  SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
+  SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst &
+#ifdef ISO
+       ,xtsnow_rst,Rland_ice_rst &
+#endif       
+       )
 
     USE indice_sol_mod
+#ifdef ISO
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    USE isotopes_verif_mod, ONLY: errmax,errmaxrel
+#endif    
+#endif
 
     INCLUDE "dimsoil.h"
@@ -3642,4 +4227,8 @@
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
     REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
+#ifdef ISO
+    REAL, DIMENSION(niso,klon, nbsrf), INTENT(OUT)     :: xtsnow_rst
+    REAL, DIMENSION(niso,klon), INTENT(OUT)            :: Rland_ice_rst
+#endif
 
  
@@ -3652,4 +4241,8 @@
     qsurf_rst(:,:)    = qsurf(:,:)
     ftsoil_rst(:,:,:) = ftsoil(:,:,:)
+#ifdef ISO
+    xtsnow_rst(:,:,:)  = xtsnow(:,:,:) 
+    Rland_ice_rst(:,:) = Rland_ice(:,:)
+#endif
 
 !****************************************************************************************
@@ -3664,4 +4257,9 @@
     IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
     IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
+#ifdef ISO
+    IF (ALLOCATED(xtsnow)) DEALLOCATE(xtsnow)
+    IF (ALLOCATED(Rland_ice)) DEALLOCATE(Rland_ice)
+    IF (ALLOCATED(Roce)) DEALLOCATE(Roce)
+#endif
 
 !jyg<
@@ -3682,5 +4280,9 @@
   SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
        evap, z0m, z0h, agesno,                                  &
-       tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)  
+       tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke &
+#ifdef ISO
+      ,xtevap  &
+#endif
+&      )  
     !albedo SB <<<
     ! Give default values where new fraction has appread
@@ -3711,4 +4313,7 @@
     REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT)        :: z0m,z0h
     REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,nbsrf), INTENT(INOUT)        :: xtevap
+#endif
 
 ! Local variables
@@ -3718,4 +4323,7 @@
     CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
     INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
+#ifdef ISO
+    INTEGER           :: ixt
+#endif
 !
 ! All at once !! 
@@ -3763,4 +4371,9 @@
                 u10m(i,nsrf)  = u10m(i,nsrf_comp1)
                 v10m(i,nsrf)  = v10m(i,nsrf_comp1)
+#ifdef ISO
+                DO ixt=1,ntraciso
+                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp1)       
+                ENDDO       
+#endif
                 IF (iflag_pbl > 1) THEN
                  tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
@@ -3818,4 +4431,10 @@
                 u10m(i,nsrf)  = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
                 v10m(i,nsrf)  = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
+#ifdef ISO
+                DO ixt=1,ntraciso
+                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) &
+                                     + xtevap(ixt,i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) 
+                ENDDO       
+#endif
                 IF (iflag_pbl > 1) THEN
                  tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
@@ -3830,4 +4449,7 @@
              agesno(i,nsrf)   = 0.
              ftsoil(i,:,nsrf) = tsurf(i,nsrf)
+#ifdef ISO            
+             xtsnow(:,i,nsrf) = 0.
+#endif
           ELSE
              pfois(nsrf) = pfois(nsrf)+ 1
