Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_parallel/xios_orchidee.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_parallel/xios_orchidee.f90	(revision 7709)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_parallel/xios_orchidee.f90	(revision 7710)
@@ -302,5 +302,4 @@
           
        ELSE IF (grid_type==regular_xy ) THEN
-
           ! Global domain
           CALL xios_set_domain_attr("domain_landpoints", ni_glo=iim_g, nj_glo=jjm_g)
@@ -414,4 +413,9 @@
           CALL xios_set_field_attr("netirrig",enabled=.FALSE.)
           CALL xios_set_field_attr("SurfStor",enabled=.FALSE.)
+          CALL xios_set_field_attr("htutempmon",enabled=.FALSE.)
+          CALL xios_set_field_attr("streamlimit",enabled=.FALSE.)
+          CALL xios_set_field_attr("StreamT_TotTend",enabled=.FALSE.)
+          CALL xios_set_field_attr("StreamT_AdvTend",enabled=.FALSE.)
+          CALL xios_set_field_attr("StreamT_RelTend",enabled=.FALSE.)
        END IF
 
@@ -433,4 +437,5 @@
           CALL xios_set_field_attr("floodmap",enabled=.FALSE.)
           CALL xios_set_field_attr("floodh",enabled=.FALSE.)       
+          CALL xios_set_field_attr("floodr",enabled=.FALSE.)
           CALL xios_set_field_attr("floodout",enabled=.FALSE.)
           CALL xios_set_field_attr("flood_frac",enabled=.FALSE.)       
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_highres.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_highres.f90	(revision 7709)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_highres.f90	(revision 7710)
@@ -1,3 +1,3 @@
-! =================================================================================================================================
+! ================================================================================================================================
 ! MODULE       : routing_highres
 !
@@ -96,4 +96,9 @@
   INTEGER(i_std), SAVE                                       :: overflow_repetition = 1     !! Number of repetition of overflow for each routing step
   !$OMP THREADPRIVATE(overflow_repetition)
+  ! Soil temperature depth to be used to estimate runoff and drainage temperatures
+  !
+  REAL(r_std), PARAMETER, DIMENSION(2) :: runofftempdepth = (/ 0.0, 0.3 /)                  !! Layer which will determine the temperature of runoff
+  REAL(r_std), PARAMETER, DIMENSION(2) :: drainagetempdepth = (/ 3.0, 90.0 /)               !! Layer which will determine the temperature of runoff
+  !
   !
   !  Relation between maximum surface of ponds and basin surface, and drainage (mm/j) to the slow_res
@@ -123,6 +128,4 @@
   LOGICAL, SAVE                                              :: doponds = .FALSE.           !! Logical to choose if ponds are activated or not (true/false)
 !$OMP THREADPRIVATE(doponds)
-  LOGICAL, SAVE                                              :: do_rivertemp = .FALSE.
-!$OMP THREADPRIVATE(do_rivertemp)
   REAL(r_std), SAVE                                          :: conduct_factor = 1.         !! Adjustment factor for floodplains reinfiltration
 !$OMP THREADPRIVATE(conduct_factor)
@@ -221,6 +224,6 @@
 !!!
 
-  REAL(r_std),  SAVE, ALLOCATABLE, DIMENSION(:,:)            :: stempdiag_mean              !! Averaged soil temperatures
-!$OMP THREADPRIVATE(stempdiag_mean)
+  REAL(r_std),  SAVE, ALLOCATABLE, DIMENSION(:,:)            :: tempdiag_mean              !! Averaged soil temperatures
+!$OMP THREADPRIVATE(tempdiag_mean)
 !
 ! FLOOD OVERFLOW
@@ -404,5 +407,5 @@
   SUBROUTINE routing_highres_initialize( kjit,       nbpt,           index,                 &
                                 rest_id,     hist_id,        hist2_id,   lalo,      &
-                                neighbours,  resolution,     contfrac,   stempdiag, &
+                                neighbours,  resolution,     contfrac,   tempdiag, &
                                 returnflow,  reinfiltration, irrigation, riverflow, &
                                 coastalflow, flood_frac,     flood_res )
@@ -423,5 +426,5 @@
     REAL(r_std), INTENT(in)        :: resolution(nbpt,2)   !! The size of each grid box in X and Y (m)
     REAL(r_std), INTENT(in)        :: contfrac(nbpt)       !! Fraction of land in each grid box (unitless;0-1)
-    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
+    REAL(r_std), INTENT(in)        :: tempdiag(nbpt,ngrnd) !! Diagnostic soil temperature profile
 
     !! 0.2 Output variables
@@ -455,5 +458,5 @@
 
     CALL routing_hr_init (kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
-         riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id)
+         riverflow, coastalflow, flood_frac, flood_res, tempdiag, rest_id)
 
     routing_area => routing_area_loc  
@@ -643,8 +646,4 @@
     IF (printlev >= 5) WRITE(numout,*) 'End of routing_highres_initialize'
 
-    !! Define XIOS axis size needed for the model output
-    !CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/))
-    !CALL xios_orchidee_addaxis("nbasmon", nbasmon, (/(REAL(ib,r_std),ib=1,nbasmon)/))
-    
   END SUBROUTINE routing_highres_initialize
 
@@ -758,5 +757,5 @@
        & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
        & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
-       & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
+       & tempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
 
     IMPLICIT NONE
@@ -782,5 +781,5 @@
     REAL(r_std), INTENT(in)        :: k_litt(nbpt)         !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
     REAL(r_std), INTENT(in)        :: humrel(nbpt,nvm)     !! Soil moisture stress, root extraction potential (unitless)
-    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
+    REAL(r_std), INTENT(in)        :: tempdiag(nbpt,ngrnd) !! Diagnostic soil temperature profile
     REAL(r_std), INTENT(in)        :: reinf_slope(nbpt)    !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
 
@@ -816,6 +815,4 @@
     REAL(r_std), DIMENSION(nbpt,nbasmax)   :: stemp_total_tend, stemp_advec_tend, stemp_relax_tend
     !
-    LOGICAL, SAVE                  :: xios_sendonce=.TRUE.
-    !
 !_ ================================================================================================================================
 
@@ -834,5 +831,5 @@
     runoff_mean(:) = runoff_mean(:) + runoff(:)
     drainage_mean(:) = drainage_mean(:) + drainage(:)
-    floodtemp(:) = stempdiag(:,floodtemp_lev)
+    floodtemp(:) = tempdiag(:,floodtemp_lev)
     precip_mean(:) =  precip_mean(:) + precip_rain(:)
     !
@@ -867,5 +864,5 @@
     totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dt_sechiba/dt_routing
     k_litt_mean(:) = k_litt_mean(:) + k_litt(:)*dt_sechiba/dt_routing
-    stempdiag_mean(:,:) = stempdiag_mean(:,:) + stempdiag(:,:)*dt_sechiba/dt_routing
+    tempdiag_mean(:,:) = tempdiag_mean(:,:) + tempdiag(:,:)*dt_sechiba/dt_routing
     !
     ! Only potentially vegetated surfaces are taken into account. At the start of
@@ -889,5 +886,5 @@
        CALL routing_hr_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, &
             & vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, &
-            & stempdiag_mean, reinf_slope, lakeinflow_mean, returnflow_mean, reinfiltration_mean, &
+            & tempdiag_mean, reinf_slope, lakeinflow_mean, returnflow_mean, reinfiltration_mean, &
             & irrigation_mean, riverflow_mean, coastalflow_mean, hydrographs, slowflow_diag, flood_frac, &
             & flood_res, netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, &
@@ -911,5 +908,5 @@
        totnobio_mean(:) = zero
        k_litt_mean(:) = zero
-       stempdiag_mean(:,:) = zero
+       tempdiag_mean(:,:) = zero
        vegtot_mean(:) = zero
 
@@ -937,9 +934,7 @@
        CALL xios_orchidee_send_field("wbr_lake",  (lake_diag   - lake_diag_old - &
             lakeinflow_mean + return_lakes)/dt_routing)
-       IF ( do_rivertemp ) THEN
-          CALL xios_orchidee_send_field("StreamT_TotTend", stemp_total_tend)
-          CALL xios_orchidee_send_field("StreamT_AdvTend", stemp_advec_tend)
-          CALL xios_orchidee_send_field("StreamT_RelTend", stemp_relax_tend)
-       ENDIF
+       CALL xios_orchidee_send_field("StreamT_TotTend", stemp_total_tend)
+       CALL xios_orchidee_send_field("StreamT_AdvTend", stemp_advec_tend)
+       CALL xios_orchidee_send_field("StreamT_RelTend", stemp_relax_tend)
     ENDIF
 
@@ -956,27 +951,24 @@
     ! Write diagnostics
     !
-    IF ( xios_sendonce ) THEN
-       !
-       CALL xios_orchidee_send_field("mask_coast",mask_coast)
-
-       IF ( do_irrigation ) THEN 
-          CALL xios_orchidee_send_field("irrigmap",irrigated)
-       ENDIF
+    !
+    CALL xios_orchidee_send_field("mask_coast",mask_coast)
+
+    IF ( do_irrigation ) THEN 
+       CALL xios_orchidee_send_field("irrigmap",irrigated)
+    ENDIF
        
-       IF ( do_floodplains ) THEN
-          !! May be improved by performing the operation with XIOS
-          floodmap(:) = 0.0
-          DO ig=1,nbpt
-             floodmap(ig) = SUM(floodplains(ig,:)) / (area(ig)*contfrac(ig))
-          END DO
-          CALL xios_orchidee_send_field("floodmap",floodmap)
-       ENDIF
+    IF ( do_floodplains ) THEN
+       !! May be improved by performing the operation with XIOS
+       floodmap(:) = 0.0
+       DO ig=1,nbpt
+          floodmap(ig) = SUM(floodplains(ig,:)) / (area(ig)*contfrac(ig))
+       END DO
+       CALL xios_orchidee_send_field("floodmap",floodmap)
+    ENDIF
        
-       IF ( doswamps ) THEN
-          CALL xios_orchidee_send_field("swampmap",swamp)
-       ENDIF
+    IF ( doswamps ) THEN
+       CALL xios_orchidee_send_field("swampmap",swamp)
+    ENDIF
        
-       xios_sendonce=.FALSE.
-    ENDIF
     !
     ! Water storage in reservoirs [kg/m^2]
@@ -1009,9 +1001,8 @@
     CALL xios_orchidee_send_field("hydrographs",hydrographs/mille/dt_sechiba)
     CALL xios_orchidee_send_field("htuhgmon",HTUhgmon/mille/dt_sechiba)
-    IF ( do_rivertemp ) THEN
-       CALL xios_orchidee_send_field("htutempmon",HTUtempmon)
-       CALL xios_orchidee_send_field("hydrotemp", hydrotemp)
-       CALL xios_orchidee_send_field("streamlimit", streamlimit)
-    ENDIF
+    CALL xios_orchidee_send_field("htutempmon",HTUtempmon)
+    CALL xios_orchidee_send_field("hydrotemp", hydrotemp)
+    CALL xios_orchidee_send_field("streamlimit", streamlimit)
+
     CALL xios_orchidee_send_field("slowflow",slowflow_diag/mille/dt_sechiba) ! previous id name: Qb
     CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
@@ -1164,9 +1155,7 @@
     CALL restput_p (rest_id, 'flood_res', nbp_glo, 1, 1, kjit, flood_res, 'scatter', nbp_glo, index_g)
 
-    IF ( do_rivertemp ) THEN
-       CALL restput_p (rest_id, 'fasttemp', nbp_glo, nbasmax, 1, kjit, fast_temp, 'scatter',  nbp_glo, index_g)
-       CALL restput_p (rest_id, 'slowtemp', nbp_glo, nbasmax, 1, kjit, slow_temp, 'scatter',  nbp_glo, index_g)
-       CALL restput_p (rest_id, 'streamtemp', nbp_glo, nbasmax, 1, kjit, stream_temp, 'scatter',nbp_glo,index_g)
-    ENDIF
+    CALL restput_p (rest_id, 'fasttemp', nbp_glo, nbasmax, 1, kjit, fast_temp, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'slowtemp', nbp_glo, nbasmax, 1, kjit, slow_temp, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'streamtemp', nbp_glo, nbasmax, 1, kjit, stream_temp, 'scatter',nbp_glo,index_g)
  
     
@@ -1182,8 +1171,6 @@
     CALL restput_p (rest_id, 'htuhgmon', nbp_glo, nbasmon, 1, kjit, HTUhgmon, 'scatter',  nbp_glo, index_g)
     CALL restput_p (rest_id, 'slowflow_diag', nbp_glo, 1, 1, kjit, slowflow_diag, 'scatter',  nbp_glo, index_g)
-    IF ( do_rivertemp ) THEN
-       CALL restput_p (rest_id, 'hydrotemp', nbp_glo, 1, 1, kjit, hydrotemp, 'scatter',  nbp_glo, index_g)
-       CALL restput_p (rest_id, 'htutempmon', nbp_glo, nbasmon, 1, kjit, HTUtempmon, 'scatter',  nbp_glo, index_g)
-    ENDIF
+    CALL restput_p (rest_id, 'hydrotemp', nbp_glo, 1, 1, kjit, hydrotemp, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'htutempmon', nbp_glo, nbasmon, 1, kjit, HTUtempmon, 'scatter',  nbp_glo, index_g)
     !
     ! Keep track of the accumulated variables
@@ -1198,5 +1185,5 @@
     CALL restput_p (rest_id, 'k_litt_route', nbp_glo, 1, 1, kjit, k_litt_mean, 'scatter',  nbp_glo, index_g)
     CALL restput_p (rest_id, 'vegtot_route', nbp_glo, 1, 1, kjit, vegtot_mean, 'scatter',  nbp_glo, index_g)
-    CALL restput_p (rest_id, 'stempdiag_route', nbp_glo, nslm, 1, kjit, stempdiag_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'tempdiag_route', nbp_glo, ngrnd, 1, kjit, tempdiag_mean, 'scatter',  nbp_glo, index_g)
 
     CALL restput_p (rest_id, 'gridrephtu', nbp_glo, 1, 1, kjit, REAL(hydrodiag,r_std), 'scatter',  nbp_glo, index_g)
@@ -1240,5 +1227,5 @@
 
   SUBROUTINE routing_hr_init(kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
-       &                  riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id)
+       &                  riverflow, coastalflow, flood_frac, flood_res, tempdiag, rest_id)
     !
     IMPLICIT NONE
@@ -1250,5 +1237,5 @@
     INTEGER(i_std), INTENT(in)                   :: nbpt           !! Domain size (unitless)
     INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index          !! Indices of the points on the map (unitless)
-    REAL(r_std), DIMENSION(nbpt,nslm),INTENT(in) :: stempdiag      !! Temperature profile in soil
+    REAL(r_std), DIMENSION(nbpt,ngrnd),INTENT(in) :: tempdiag      !! Temperature profile in soil
     INTEGER(i_std), INTENT(in)                   :: rest_id        !! Restart file identifier (unitless)
     !
@@ -1287,5 +1274,5 @@
     !Config Units = [seconds]
     !
-    dt_routing = one_day
+    dt_routing = dt_sechiba
     CALL getin_p('DT_ROUTING', dt_routing)
     !
@@ -1442,12 +1429,4 @@
     CALL getin_p("MAX_LAKE_RESERVOIR", max_lake_reservoir)
 
-    !Config Key   = DO_RIVERTEMP
-    !Config Desc  = Activates the river temperature calculations
-    !Config If    = 
-    !Config Def   = False
-    !Config Help  = 
-    !Config Units = -
-    do_rivertemp = .FALSE.
-    CALL getin_p("DO_RIVERTEMP", do_rivertemp)
     !
     !
@@ -1963,35 +1942,33 @@
     irrigation(:) = irrigation_mean(:)
     
-    IF ( do_rivertemp ) THEN
-       ALLOCATE (fast_temp(nbpt,nbasmax), stat=ier)
-       IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fast_temp','','')
-       var_name = 'fasttemp'
-       CALL ioconf_setatt_p('UNITS', 'K')
-       CALL ioconf_setatt_p('LONG_NAME','Water temperature in the fast reservoir')
-       CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_temp, "gather", nbp_glo, index_g)
-
-       ALLOCATE (slow_temp(nbpt,nbasmax), stat=ier)
-       IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for slow_temp','','')
-       var_name = 'slowtemp'
-       CALL ioconf_setatt_p('UNITS', 'K')
-       CALL ioconf_setatt_p('LONG_NAME','Water temperature in the slow reservoir')
-       CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_temp, "gather", nbp_glo, index_g)
+    ALLOCATE (fast_temp(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fast_temp','','')
+    var_name = 'fasttemp'
+    CALL ioconf_setatt_p('UNITS', 'K')
+    CALL ioconf_setatt_p('LONG_NAME','Water temperature in the fast reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_temp, "gather", nbp_glo, index_g)
+
+    ALLOCATE (slow_temp(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for slow_temp','','')
+    var_name = 'slowtemp'
+    CALL ioconf_setatt_p('UNITS', 'K')
+    CALL ioconf_setatt_p('LONG_NAME','Water temperature in the slow reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_temp, "gather", nbp_glo, index_g)
        
-       IF ( COUNT(fast_temp == val_exp) == nbpt*nbasmax ) THEN
-          CALL groudwatertemp(nbpt, nbasmax, nslm, stempdiag, diaglev, fast_temp, slow_temp)
-       ENDIF
+    IF ( COUNT(fast_temp == val_exp) == nbpt*nbasmax ) THEN
+       CALL groundwatertemp(nbpt, nbasmax, ngrnd, tempdiag, znt, dlt, fast_temp, slow_temp)
+    ENDIF
        
-       ALLOCATE (stream_temp(nbpt,nbasmax), stat=ier)
-       IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_temp','','')
-       var_name = 'streamtemp'
-       CALL ioconf_setatt_p('UNITS', 'K')
-       CALL ioconf_setatt_p('LONG_NAME','Water temperature in the stream reservoir')
-       CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_temp, "gather", nbp_glo, index_g)
-
-       IF ( COUNT(stream_temp == val_exp) == nbpt*nbasmax ) THEN
-          DO ig=1,nbpt 
-             stream_temp(ig,:) = stempdiag(ig,1)
-          ENDDO
-       ENDIF
+    ALLOCATE (stream_temp(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_temp','','')
+    var_name = 'streamtemp'
+    CALL ioconf_setatt_p('UNITS', 'K')
+    CALL ioconf_setatt_p('LONG_NAME','Water temperature in the stream reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_temp, "gather", nbp_glo, index_g)
+
+    IF ( COUNT(stream_temp == val_exp) == nbpt*nbasmax ) THEN
+       DO ig=1,nbpt 
+         stream_temp(ig,:) = tempdiag(ig,1)
+       ENDDO
     ENDIF
     
@@ -2019,5 +1996,5 @@
     ALLOCATE (floodtemp(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodtemp','','')
-    floodtemp(:) = stempdiag(:,floodtemp_lev)
+    floodtemp(:) = tempdiag(:,floodtemp_lev)
     
     ALLOCATE(hydrographs(nbpt), stat=ier)
@@ -2244,11 +2221,11 @@
     CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un)
     !
-    ALLOCATE(stempdiag_mean(nbpt,nslm), stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stempdiag_mean','','')
-    var_name = 'stempdiag_route'
+    ALLOCATE(tempdiag_mean(nbpt,ngrnd), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for tempdiag_mean','','')
+    var_name = 'tempdiag_route'
     CALL ioconf_setatt_p('UNITS', 'K')
     CALL ioconf_setatt_p('LONG_NAME','Mean temperature profile')
-    CALL restget_p (rest_id, var_name, nbp_glo, nslm, 1, kjit, .TRUE., stempdiag_mean, "gather", nbp_glo, index_g)
-    CALL setvar_p (stempdiag_mean, val_exp, 'NO_KEYWORD', Zero)
+    CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., tempdiag_mean, "gather", nbp_glo, index_g)
+    CALL setvar_p (tempdiag_mean, val_exp, 'NO_KEYWORD', Zero)
     !
     DEALLOCATE(tmp_real_g)
@@ -2355,5 +2332,5 @@
     IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_mean)
     IF (ALLOCATED(k_litt_mean)) DEALLOCATE(k_litt_mean)
-    IF (ALLOCATED(stempdiag_mean)) DEALLOCATE(stempdiag_mean)
+    IF (ALLOCATED(tempdiag_mean)) DEALLOCATE(tempdiag_mean)
     IF (ALLOCATED(totnobio_mean)) DEALLOCATE(totnobio_mean)
     IF (ALLOCATED(vegtot_mean)) DEALLOCATE(vegtot_mean)
@@ -2466,5 +2443,5 @@
 
   SUBROUTINE routing_hr_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, &
-       &                  vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, stempdiag, &
+       &                  vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, tempdiag, &
        &                  reinf_slope, lakeinflow, returnflow, reinfiltration, irrigation, riverflow, &
        &                  coastalflow, hydrographs, slowflow_diag, flood_frac, flood_res, &
@@ -2488,5 +2465,5 @@
     REAL(r_std), INTENT(in)                      :: k_litt(nbpt)              !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
     REAL(r_std), INTENT(in)                      :: floodtemp(nbpt)           !! Temperature to decide if floodplains work (K)
-    REAL(r_std), INTENT(in)                      :: stempdiag(nbpt,nslm)      !! Soil temperature profiles (K)
+    REAL(r_std), INTENT(in)                      :: tempdiag(nbpt,ngrnd)      !! Soil temperature profiles (K)
     REAL(r_std), INTENT(in)                      :: reinf_slope(nbpt)         !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
     REAL(r_std), INTENT(out)                     :: lakeinflow(nbpt)          !! Water inflow to the lakes (kg/dt)
@@ -2521,5 +2498,5 @@
     REAL(r_std), DIMENSION(nbpt, nbasmax)        :: return_swamp              !! Inflow to the swamp (kg/dt)
     REAL(r_std), DIMENSION(nbpt, nbasmax)        :: source
-    REAL(r_std), DIMENSION(nbpt, nbasmax)        :: ewh, wrr
+    REAL(r_std), DIMENSION(nbpt, nbasmax)        :: ewh
     !
     ! Irrigation per basin
@@ -2580,26 +2557,19 @@
 
     REAL(r_std)                                  :: reduced                   !! Discharge reduction due to floodplains
-    REAL(r_std)                                  :: restime, minrestime, maxrestime!! Scaling for residence times in the relaxation of temperature.
     REAL(r_std)                                  :: htmp, hscale              !! Water height scalingfor temperature relaxation
     REAL(r_std)                                  :: krelax, den
     !! PARAMETERS
     LOGICAL, PARAMETER                           :: check_reservoir = .FALSE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false)
-    LOGICAL, SAVE                                :: origformula = .FALSE.
 !_ ================================================================================================================================
     !
-    minrestime = 86400/2.0
-    maxrestime = stream_maxresid * stream_tcst
-    !
-    origformula = .FALSE.
-    CALL getin_p('ELIOTTFORMULA',origformula)
+    !
     hscale = 1.
-    CALL getin_p('HSCALEKH',hscale)
+    CALL getin_p('ROUTING_HSCALEKH',hscale)
     !
     transport(:,:) = zero
     transport_glo(:,:) = zero
-    IF ( do_rivertemp ) THEN
-       transport_temp(:,:) = zero !tp_00  its a transport, not a temperature !!
-       transport_temp_glo(:,:) = zero !tp_00
-    ENDIF
+    transport_temp(:,:) = zero !tp_00  its a transport, not a temperature !!
+    transport_temp_glo(:,:) = zero !tp_00
+    
     irrig_netereq(:) = zero
     irrig_needs(:,:) = zero
@@ -2631,7 +2601,5 @@
 !> the sub-basin considered and its downstream neighbor.
     !
-    IF ( do_rivertemp ) THEN
-       CALL groudwatertemp(nbpt, nbasmax, nslm, stempdiag, diaglev, fast_temp, slow_temp)
-    ENDIF
+    CALL groundwatertemp(nbpt, nbasmax, ngrnd, tempdiag, znt, dlt, fast_temp, slow_temp)
     !
     streamlimit(:) = zero
@@ -2799,9 +2767,7 @@
     source(:,:) = fast_flow(:,:) + slow_flow(:,:) + stream_flow(:,:)
     CALL downstreamsum(nbpt, nbasmax, source, transport)
-    IF ( do_rivertemp ) THEN
-       source(:,:) = fast_flow(:,:)*fast_temp(:,:) + slow_flow(:,:)*slow_temp(:,:) +  &
+    source(:,:) = fast_flow(:,:)*fast_temp(:,:) + slow_flow(:,:)*slow_temp(:,:) +  &
             &                  stream_flow(:,:)*stream_temp(:,:)
-       CALL downstreamsum(nbpt, nbasmax, source, transport_temp)
-    ENDIF
+    CALL downstreamsum(nbpt, nbasmax, source, transport_temp)
     !-
     !- Do the floodings - First initialize
@@ -2868,69 +2834,41 @@
                & slow_flow(ig,ib)
           !
-          IF ( do_rivertemp ) THEN
-             oldstream = stream_reservoir(ig, ib) * stream_temp(ig,ib)
-          ENDIF
+          oldstream = stream_reservoir(ig, ib) * stream_temp(ig,ib)
           !
           stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_flow(ig,ib) + transport(ig,ib) - &
                & stream_flow(ig,ib) - return_swamp(ig,ib) - floods(ig,ib)
           !
-          ! Code d'Eliott
-          !
-          IF ( do_rivertemp ) THEN
+          ! Diagnostics of the stream reservoir 
+          !
+          IF ( routing_area(ig,ib) > zero ) THEN
+             ! 1000 to transform kg into m^3
+             htmp = stream_reservoir(ig,ib)*1000/routing_area(ig,ib)
+             ewh(ig,ib) = 1.0/(1.0+htmp*hscale)
+          ELSE
+             ewh(ig,ib) = 1.0
+          ENDIF
+          !
+          !reste du calcul
+          !
+          krelax = ewh(ig,ib)
+          !
+          den = 1.0/(1.0+dt_routing*krelax)
+          IF ( stream_reservoir(ig,ib) > 1.e-6 ) THEN
+             oldtemp = stream_temp(ig,ib)
+             stream_temp(ig,ib) = den * dt_routing * krelax * fast_temp(ig,ib) + &
+                  & den * oldstream/stream_reservoir(ig,ib) + &
+                  & den * transport_temp(ig, ib)/stream_reservoir(ig,ib) - &
+                  & den * oldtemp*stream_flow(ig,ib)/stream_reservoir(ig,ib)
              !
-             ! Diagnostics of the stream reservoir 
-             !
-             IF ( routing_area(ig,ib) > zero ) THEN
-                ! 1000 to transform kg into m^3
-                htmp = stream_reservoir(ig,ib)*1000/routing_area(ig,ib)
-                ewh(ig,ib) = 1.0/(1.0+htmp*hscale)
-             ELSE
-                ewh(ig,ib) = 1.0
-             ENDIF
-             !
-             IF (stream_reservoir(ig,ib) > zero ) THEN
-                ! Residence time in seconds
-                restime = stream_resid(ig,ib)*stream_tcst
-                wrr(ig,ib) =ABS((restime-minrestime)/(maxrestime-minrestime))
-             ELSE
-                wrr(ig,ib) = 1.0
-             ENDIF
-             !
-             !reste du calcul
-             !
-             IF ( origformula ) THEN
-                IF (stream_reservoir(ig,ib) .LT. 1.) THEN
-                   !
-                   stream_temp(ig,ib) = max(stempdiag(ig,1), ZeroCelsius)
-                   ! 
-                ELSE
-                   !
-                   stream_temp(ig,ib) = oldstream/stream_reservoir(ig,ib) + &
-                        & transport_temp(ig, ib)/stream_reservoir(ig,ib) - &
-                        & stream_temp(ig,ib)*stream_flow(ig,ib)/stream_reservoir(ig,ib)
-                ENDIF
-             ELSE
-                krelax = ewh(ig,ib)
-                !
-                den = 1.0/(1.0+dt_routing*krelax)
-                IF ( stream_reservoir(ig,ib) > 1.e-6 ) THEN
-                   oldtemp = stream_temp(ig,ib)
-                   stream_temp(ig,ib) = den * dt_routing * krelax * fast_temp(ig,ib) + &
-                        & den * oldstream/stream_reservoir(ig,ib) + &
-                        & den * transport_temp(ig, ib)/stream_reservoir(ig,ib) - &
-                        & den * oldtemp*stream_flow(ig,ib)/stream_reservoir(ig,ib)
-                   !
-                   !Stream_temp [K], stream_reservoir [kg], WaterCp [J/g/K] yields tendencies in GJ/s
-                   ! 
-                   stemp_total_tend(ig,ib) = WaterCp*1.e-6*(stream_temp(ig,ib)*stream_reservoir(ig,ib) - oldstream)/dt_routing
-                   stemp_advec_tend(ig,ib) = WaterCp*1.e-6*(transport_temp(ig, ib) - oldtemp*stream_flow(ig,ib))/dt_routing
-                   stemp_relax_tend(ig,ib) = WaterCp*1.e-6*stream_reservoir(ig,ib)*krelax*(fast_temp(ig,ib)-stream_temp(ig,ib))
-                ELSE
-                   stream_temp(ig,ib) = MAX(fast_temp(ig,ib), ZeroCelsius)
-                   stemp_total_tend(ig,ib) = zero
-                   stemp_advec_tend(ig,ib) = zero
-                   stemp_relax_tend(ig,ib) = zero
-                ENDIF
-             ENDIF
+             !Stream_temp [K], stream_reservoir [kg], WaterCp [J/g/K] yields tendencies in GJ/s
+             ! 
+             stemp_total_tend(ig,ib) = WaterCp*1.e-6*(stream_temp(ig,ib)*stream_reservoir(ig,ib) - oldstream)/dt_routing
+             stemp_advec_tend(ig,ib) = WaterCp*1.e-6*(transport_temp(ig, ib) - oldtemp*stream_flow(ig,ib))/dt_routing
+             stemp_relax_tend(ig,ib) = WaterCp*1.e-6*stream_reservoir(ig,ib)*krelax*(fast_temp(ig,ib)-stream_temp(ig,ib))
+          ELSE
+             stream_temp(ig,ib) = MAX(fast_temp(ig,ib), ZeroCelsius)
+             stemp_total_tend(ig,ib) = zero
+             stemp_advec_tend(ig,ib) = zero
+             stemp_relax_tend(ig,ib) = zero
           ENDIF
           !
@@ -3264,7 +3202,5 @@
              IF (HTUdiag_loc(ig,im) > 0 .AND. HTUdiag_loc(ig,im) .EQ. ib ) THEN
                 HTUhgmon(ig,im) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib)
-                IF ( do_rivertemp ) THEN
-                   HTUtempmon(ig,im) = stream_temp(ig,ib)
-                ENDIF
+                HTUtempmon(ig,im) = stream_temp(ig,ib)
              ENDIF
           ENDDO
@@ -3272,7 +3208,5 @@
           IF (hydrodiag(ig) == ib) THEN
              hydrographs(ig) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib)
-             IF ( do_rivertemp ) THEN
-                hydrotemp(ig) = stream_temp(ig,ib)
-             ENDIF
+             hydrotemp(ig) = stream_temp(ig,ib)
              slowflow_diag(ig) = slowflow_diag(ig) + slow_flow(ig,ib)
           ENDIF
@@ -3337,5 +3271,5 @@
   !
 !! ================================================================================================================================
-!! SUBROUTINE 	: groundwater_temp
+!! SUBROUTINE 	: groundwatertemp
 !!
 !>\BRIEF        : This subroutine computes the temperature of the groundwater leaving the HTU
@@ -3354,54 +3288,84 @@
 !_ ================================================================================================================================
     !-
-  SUBROUTINE groudwatertemp(nbpt, nbasmax, nslm, stempdiag, diaglev, fast_temp, slow_temp)
+  SUBROUTINE groundwatertemp(nbpt, nbasmax, nl, tempdiag, lev, dlz, fast_temp, slow_temp)
     ! INPUT
-    INTEGER(i_std), INTENT(in)                   :: nbpt, nbasmax, nslm
-    REAL(r_std), INTENT(in)                      :: stempdiag(nbpt,nslm)
-    REAL(r_std), INTENT(in)                      :: diaglev(nslm)
+    INTEGER(i_std), INTENT(in)                   :: nbpt, nbasmax, nl
+    REAL(r_std), INTENT(in)                      :: tempdiag(nbpt,nl)
+    REAL(r_std), INTENT(in)                      :: lev(nl), dlz(nl)
     REAL(r_std), INTENT(inout)                   :: slow_temp(nbpt,nbasmax), fast_temp(nbpt,nbasmax)
     ! OUTPUT
     ! LOCAL
     INTEGER(i_std)                               :: ig, ib, im
-    REAL(r_std)                                  :: ft, st, fw, sw, ubnd, lbnd, dz
-    LOGICAL, SAVE                                :: lowestt = .FALSE., alltop=.FALSE.
-    !
-    lowestt=.FALSE.
-    CALL getin_p('LOWESTT', lowestt)
-    alltop=.FALSE.
-    CALL getin_p('ALLTOPTT', alltop)
-    !
-    DO ib=1,nbasmax
-       DO ig=1,nbpt
-          ft=zero
-          st=zero
-          fw=zero
-          sw=zero
-          ubnd=zero
-          lbnd=diaglev(nslm)+(diaglev(nslm)+diaglev(nslm-1))*0.5
-          ! The sums are weighted by the thickness of the soil layers !
-          DO im=1,ntemp_layer
-             dz = ((diaglev(im)+diaglev(im+1))*0.5-ubnd)
-             ft = ft + stempdiag(ig,im)*dz
-             fw = fw + dz
-             ubnd=(diaglev(im)+diaglev(im+1))*0.5
-             !
-             dz = (lbnd-(diaglev(nslm-(im))+diaglev(nslm-(im)+1))*0.5)
-             st = st + stempdiag(ig,nslm-(im-1))*dz
-             sw = sw + dz
-             lbnd = (diaglev(nslm-(im))+diaglev(nslm-(im)+1))*0.5
+    REAL(r_std)                                  :: sw
+    REAL(r_std)                                  :: rw(nl), dw(nl)
+    LOGICAL, SAVE                                :: alltop=.FALSE.
+    LOGICAL, SAVE                                :: FirstCall=.TRUE.
+    !
+    IF ( FirstCall ) THEN
+       !Config Key   = ROUTING_ALLTOPT
+       !Config Desc  = Should drainage have the temperature of the top soil (0.3m) ?
+       !Config Def   = False
+       !Config Help  = The default behaviour of the scheme is that runoff has the temperature
+       !Config Help    of the top 30 cm of soil. Drainage will have the temperature of the lowest
+       !Config Help    soil layer (3-17m). If set to True this flag will give drainage the same
+       !Config Help    temperature as runoff.
+       !Config Units = Logical
+       alltop=.FALSE.
+       CALL getin_p('ROUTING_ALLTOPT', alltop)
+       !
+       WRITE(numout,*) "Runoff will have the average soil temperature of layers from ", runofftempdepth(1),&
+            &          " to ", runofftempdepth(2), "[m]"
+       !
+       IF ( alltop ) THEN
+          WRITE(numout,*) "Drainage will have the average soil temperature of layers from ", runofftempdepth(1),&
+            &          " to ", runofftempdepth(2), "[m]"
+       ELSE
+          WRITE(numout,*) "Drainage will have the average soil temperature of layers from ", drainagetempdepth(1),&
+               &          " to ", MIN(drainagetempdepth(2), SUM(dlz)), "[m]"
+       ENDIF
+       FirstCall=.FALSE.
+    ENDIF
+    !
+    CALL tempdepthweight(nl, dlz, runofftempdepth(1), runofftempdepth(2), rw)
+    CALL tempdepthweight(nl, dlz, drainagetempdepth(1), MIN(drainagetempdepth(2), SUM(dlz)), dw)
+    !
+    slow_temp(:,:) = zero
+    fast_temp(:,:) = zero
+    ! Compute for each HTU the temperature of runoff and drainage water.
+    DO im = 1,nl
+       DO ib=1,nbasmax
+          DO ig=1,nbpt
+             fast_temp(ig,ib) = fast_temp(ig,ib) + tempdiag(ig,im)*rw(im)
+             ! The option to have drainage water at the same temperature as runoff
+             IF ( alltop ) THEN
+                slow_temp(ig,ib) = slow_temp(ig,ib) + tempdiag(ig,im)*rw(im)
+             ELSE
+                slow_temp(ig,ib) = slow_temp(ig,ib) + tempdiag(ig,im)*dw(im)
+             ENDIF
           ENDDO
-          !
-          IF ( lowestt ) THEN
-             slow_temp(ig,ib) = stempdiag(ig,nslm)
-          ELSE
-             slow_temp(ig,ib) = st/sw
-          ENDIF
-          fast_temp(ig,ib) = ft/fw
-          IF ( alltop ) THEN
-             slow_temp(ig,ib) = fast_temp(ig,ib)
-          ENDIF
        ENDDO
     ENDDO
-  END SUBROUTINE groudwatertemp
+
+  END SUBROUTINE groundwatertemp
+
+  SUBROUTINE tempdepthweight(n, dz, top, bot, w)
+    ! Input
+    INTEGER(i_std), INTENT(in)                   :: n
+    REAL(r_std), INTENT(in)                      :: dz(n)
+    REAL(r_std), INTENT(in)                      :: top, bot
+    ! Output
+    REAL(r_std), INTENT(out)                     :: w(n)
+    ! Local
+    INTEGER(i_std)                               :: i
+    REAL(r_std)                                  :: sw
+    w(:) = zero
+    sw = zero
+    DO i=1,n
+       w(i) = MAX(zero, MIN(sw+dz(i), bot) - MAX(top, sw))
+       sw = sw + dz(i)
+    ENDDO
+    w(:) = w(:)/(bot-top)
+  END SUBROUTINE tempdepthweight
+
 !! ================================================================================================================================
 !! SUBROUTINE 	: downstreamsum
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90	(revision 7709)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90	(revision 7710)
@@ -100,9 +100,8 @@
        kjit,        nbpt,           index,                 &
        rest_id,     hist_id,        hist2_id,   lalo,      &
-       neighbours,  resolution,     contfrac,   stempdiag, &
+       neighbours,  resolution,     contfrac,   stempdiag, ftempdiag, &
        soiltile,    irrig_frac,     veget_max,  irrigated_next, &    
        returnflow,  reinfiltration, irrigation, riverflow, &
        coastalflow, flood_frac,     flood_res )
-       
 
 
@@ -121,4 +120,5 @@
     REAL(r_std), INTENT(in)        :: contfrac(nbpt)       !! Fraction of land in each grid box (unitless;0-1)
     REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
+    REAL(r_std), INTENT(in)        :: ftempdiag(nbpt,ngrnd)!! Diagnostic soil temperature profile over full column
     REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless)
     REAL(r_std), INTENT(in)        :: veget_max(nbpt,nvm)  !! Maximal fraction of vegetation (unitless;0-1) !
@@ -214,8 +214,7 @@
   SUBROUTINE routing_wrapper_main(kjit, nbpt, index, &
        lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
-       drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
-       stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, &
+       drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, stempdiag, &
+       ftempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, &
        soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) 
-
 
     IMPLICIT NONE
@@ -242,4 +241,5 @@
     REAL(r_std), INTENT(in)        :: humrel(nbpt,nvm)     !! Soil moisture stress, root extraction potential (unitless)
     REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
+    REAL(r_std), INTENT(in)        :: ftempdiag(nbpt,ngrnd)!! Diagnostic soil temperature profile over full column
     REAL(r_std), INTENT(in)        :: reinf_slope(nbpt)    !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
     REAL(r_std), INTENT(in)        :: root_deficit(nbpt)   !! soil water deficit
@@ -276,5 +276,5 @@
             lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
             drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
-            stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
+            ftempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
 
     ELSE IF(routing_method=='simple') THEN 
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 7709)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 7710)
@@ -13,5 +13,5 @@
 !! processes as well.
 !!
-!!\n DESCRIPTION  : :: shumdiag, :: litterhumdiag and :: stempdiag are not 
+!!\n DESCRIPTION  : :: shumdiag, :: litterhumdiag and :: stempdiag :: ftempdiag are not 
 !! saved in the restart file because at the first time step because they 
 !! are recalculated. However, they must be saved as they are in slowproc 
@@ -148,4 +148,6 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: stempdiag      !! Temperature which controls canopy evolution (K)
 !$OMP THREADPRIVATE(stempdiag)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: ftempdiag      !! Temperature over the full soil column for river temperature (K)
+!$OMP THREADPRIVATE(ftempdiag)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintveg       !! Water on vegetation due to interception 
                                                                      !! @tex $(kg m^{-2})$ @endtex
@@ -511,5 +513,5 @@
     CALL thermosoil_initialize (kjit, kjpindex, rest_id, mcs,  &
          temp_sol_new, snow,       shumdiag_perma,        &
-         soilcap,      soilflx,    stempdiag,             &
+         soilcap,      soilflx,    stempdiag, ftempdiag,             &
          gtemp,               &
          mc_layh,  mcl_layh,   soilmoist,       njsc ,     &
@@ -525,5 +527,5 @@
             kjit,        kjpindex,       index,                 &
             rest_id,     hist_id,        hist2_id,   lalo,      &
-            neighbours,  resolution,     contfrac,   stempdiag, &
+            neighbours,  resolution,     contfrac,   stempdiag, ftempdiag, &
             soiltile,    irrig_frac,     veget_max,  irrigated_next, &
             returnflow,  reinfiltration, irrigation, riverflow, &
@@ -786,5 +788,5 @@
          index, indexgrnd, mcs, &
          temp_sol_new, snow, soilcap, soilflx, &
-         shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
+         shumdiag_perma, stempdiag, ftempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
          snowdz,snowrho,snowtemp,gtemp,pb,&
          mc_layh, mcl_layh, soilmoist, njsc,frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
@@ -797,6 +799,6 @@
        CALL routing_wrapper_main (kjit, kjpindex, index, &
             & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
-            & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
-            & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, &
+            & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, stempdiag, &
+            & ftempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id, &
             & soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw)
     ELSE
@@ -1753,4 +1755,7 @@
     IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stempdiag','','')
 
+    ALLOCATE (ftempdiag(kjpindex, ngrnd),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ftempdiag','','')
+
     ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_flux','','')
@@ -2012,4 +2017,5 @@
     IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
     IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
+    IF ( ALLOCATED (ftempdiag)) DEALLOCATE (ftempdiag)
     IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
     IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 7709)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 7710)
@@ -229,7 +229,7 @@
   !! \n
   !_ ==============================================================================================================================
-  SUBROUTINE thermosoil_initialize (kjit,          kjpindex,   rest_id,          mcs,     &
-                                    temp_sol_new,  snow,       shumdiag_perma,            &
-                                    soilcap,       soilflx,    stempdiag,                 &
+  SUBROUTINE thermosoil_initialize (kjit,          kjpindex,   rest_id,          mcs,       &
+                                    temp_sol_new,  snow,       shumdiag_perma,              &
+                                    soilcap,       soilflx,    stempdiag,        ftempdiag, &
                                     gtemp,                   &
                                     mc_layh,       mcl_layh,   tmc_layh,        njsc,     &
@@ -263,4 +263,5 @@
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)        :: soilflx          !! apparent soil heat flux considering snow and soil surface (W m-2)
     REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)   :: stempdiag        !! temperature profile on the levels in hydrol(K)
+    REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out)  :: ftempdiag        !! temperature profile on full depth for stream temperature (K)
     REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: gtemp            !! First soil layer temperature
 
@@ -519,4 +520,5 @@
     ! Send out the temperature profile on the first nslm levels(the levels treated in hydrol)
     stempdiag(:,:) = ptn(:,1:nslm)
+    ftempdiag(:,:) = ptn(:,1:ngrnd)
     
 
@@ -583,5 +585,5 @@
        index, indexgrnd, mcs, &
        temp_sol_new, snow, soilcap, soilflx, &
-       shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
+       shumdiag_perma, stempdiag, ftempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
        snowdz,snowrho,snowtemp,gtemp,pb,&
        mc_layh, mcl_layh, tmc_layh, njsc, frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
@@ -646,4 +648,5 @@
                                                                               !! see EQ3.
     REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)   :: stempdiag        !! temperature profile @tex ($K$) @endtex
+    REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out)  :: ftempdiag        !! temperature profile @tex ($K$) @endtex
     REAL(r_std),DIMENSION (kjpindex), INTENT(inout)       :: lambda_snow      !! Coefficient of the linear extrapolation of surface temperature 
     REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow       !! Integration coefficient for snow numerical scheme
@@ -768,4 +771,7 @@
     !! Initialize output arguments to be used in sechiba
     ptnlev1(:) = ptn(:,1)
+
+    !! Write the new temperature in the full in a the diagnostic variable.
+    ftempdiag(:,:) = ptn(:,1:ngrnd)
 
     !! Surface temperature is forced to zero celcius if its value is larger than melting point
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml	(revision 7709)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml	(revision 7710)
@@ -252,5 +252,11 @@
     <field id="riverflow_cpl" name="riverflow_cpl" long_name="River flow to the oceans, variable distributed over a time period before sent to ocean." unit="m^3/s"/>
     <field id="hydrographs" name="hydrographs" long_name="Hydrographs of gridbox outflow" unit="m^3/s"/>
+    <field id="hydrotemp" name="hydrotemp" long_name="River temperature at representative gridbox outflow" unit="K"/>
+    <field id="htutempmon" name="HTUstreamtemp" long_name="River temperature at monitored HTUs" unit="K" grid_ref="grid_nbasmon"/>
     <field id="htuhgmon" name="htuhgmon" long_name="Hydrographs at monitored HTUs" unit="m^3/s" grid_ref="grid_nbasmon"/>
+    <field id="streamlimit" long_name="Number of HTU with stream limiter" unit="-"/>
+    <field id="StreamT_TotTend" name="ST_TotTend" long_name="Total Stream Energy Tendency" unit="GJ/s" grid_ref="grid_nbhtu"/>
+    <field id="StreamT_AdvTend" name="ST_AdvTend" long_name="Stream Energy Advection Tendency" unit="GJ/s" grid_ref="grid_nbhtu"/>
+    <field id="StreamT_RelTend" name="ST_RelTend" long_name="Stream Energy Relaxation Tendency" unit="GJ/s" grid_ref="grid_nbhtu"/>
     <field id="fastr" name="fastr" long_name="Fast flow reservoir" unit="kg/m^2"/>
     <field id="slowr" name="slowr" long_name="Slow flow reservoir" unit="kg/m^2"/>
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml	(revision 7709)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml	(revision 7710)
@@ -74,5 +74,8 @@
     <field field_ref="riversret" level="2"/>
     <field field_ref="hydrographs" level="2"/>
-    <field field_ref="htuhgmon" name="HTUhydrographs" grid_ref="grid_nbasmon_out" level="2"/> 
+    <field field_ref="htuhgmon" name="HTUhydrographs" grid_ref="grid_nbasmon_out" level="2"/>
+    <field field_ref="hydrotemp" name="RiverTemp" level="2"/>
+    <field field_ref="htutempmon" name="HTUTemp" grid_ref="grid_nbasmon_out" level="2"/>
+    <field field_ref="streamlimit" level="2"/>
     <field field_ref="vevapnu" name="evapnu" level="0" unit="mm/d" > @vevapnu*86400  </field>
     <field field_ref="evap_bare_lim" level="2"/>
