diff --git a/src/SURFEX/modd_topodyn.F90 b/src/ARCH_SRC/surfex/dummy_topd.F90
similarity index 67%
rename from src/SURFEX/modd_topodyn.F90
rename to src/ARCH_SRC/surfex/dummy_topd.F90
index 9338f5b83aea0be923e423b01e800a39cbd74259..7d4ada880488b807fd99bd8188a169010afcddb0 100644
--- a/src/SURFEX/modd_topodyn.F90
+++ b/src/ARCH_SRC/surfex/dummy_topd.F90
@@ -1,43 +1,32 @@
-!-------------------------------------------------------------------------------
-!     ##################
-      MODULE MODD_TOPODYN
-!     ##################
+MODULE MODD_BUDGET_COUPL_ROUT
+END MODULE MODD_BUDGET_COUPL_ROUT
+!     ######################
+!     ######################
+!     ######################      
+MODULE MODD_COUPLING_TOPD
+END MODULE MODD_COUPLING_TOPD 
+!     ######################
+!     ######################
+!     ######################        
+MODULE MODD_DUMMY_EXP_PROFILE
 !
-!!****  *MODD_TOPODYN - declaration of variables used by Topodyn
-!!
-!!    PURPOSE
-!!    -------
-!
-!!
-!!**  IMPLICIT ARGUMENTS
-!!    ------------------
-!!      None 
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!     F. Habets and K. Chancibault
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original       29/09/03
-!!      BV: modifications  2006: division in two part (some variables are
-!                            now in modd_coupling_topo_n    
-!!      BV: modifications  04/2007: addition of XTOPD_STEP and NNB_TOPD_STEP
-!
-!*       0.   DECLARATIONS
-!             ------------
+END MODULE MODD_DUMMY_EXP_PROFILE
+!     ######################
+!     ######################
+!     ######################
+MODULE MODD_TOPD_PAR
 !
+REAL,    PARAMETER :: XSTEPK = 0.05 
+INTEGER, PARAMETER :: NDIM = 20 
+INTEGER, PARAMETER :: JPCAT = 10    
+END MODULE MODD_TOPD_PAR
+
+!     ######################
+!     ######################
+!     ######################
+MODULE MODD_TOPODYN
 USE MODD_TOPD_PAR, ONLY : JPCAT
-!
-IMPLICIT NONE
-!
-!-------------------------------------------------------------------------------
-! Variables specific to Topodyn
-!
- CHARACTER(LEN=15), DIMENSION(JPCAT) :: CCAT     ! base name for topographic files
+CHARACTER(LEN=15), DIMENSION(JPCAT) :: CCAT     ! base name for topographic files
 INTEGER                             :: NNCAT    ! catchments number
 !
 INTEGER                             :: NNB_TOPD_STEP   ! number of TOPODYN time steps
@@ -57,7 +46,7 @@ INTEGER, ALLOCATABLE, DIMENSION(:,:):: NLINE    ! second index of the pixel in t
 REAL, ALLOCATABLE, DIMENSION(:,:)   :: XTANB    ! pixels topographic slope (Tan(Beta))
 REAL, ALLOCATABLE, DIMENSION(:,:)   :: XSLOP    ! pixels topographic slope/length flow
 
-!Variables ā priori inutiles
+!Variables à priori inutiles
 REAL, ALLOCATABLE, DIMENSION(:,:)   :: XDAREA   ! drainage area (aire drainee)
 
 ! Variables defining the catchments
@@ -110,8 +99,38 @@ REAL, ALLOCATABLE, DIMENSION(:)   :: XRI,XRI_PREV! recharge on ISBA grid
 REAL, ALLOCATABLE, DIMENSION(:)   :: XSRFULL! reservoir of interception for
 !TOPODYN only
 REAL, ALLOCATABLE, DIMENSION(:,:) :: XDEFT! pixel deficit
-!
-!-------------------------------------------------------------------------------------
-!
+
 END MODULE MODD_TOPODYN
+!     ######################
+!     ######################
+!     ######################
+SUBROUTINE INIT_SURF_TOPD(HPROGRAM,KI)
+CHARACTER(LEN=*),  INTENT(IN)     :: HPROGRAM      !
+INTEGER,           INTENT(IN)     :: KI            ! grid dimension
+END SUBROUTINE INIT_SURF_TOPD
+
+      SUBROUTINE ISBA_TO_TOPD(PVARI,PVART)
+REAL, DIMENSION(:), INTENT(IN)      :: PVARI   ! variable from ISBA grid
+REAL, DIMENSION(:,:), INTENT(OUT)   :: PVART   ! variable for TOPODYN grid
+END SUBROUTINE ISBA_TO_TOPD
+!     ######################
+!     ######################
+!     ######################
+SUBROUTINE COUPLING_SURF_TOPD (HPROGRAM,KI)
+ CHARACTER(LEN=6), INTENT(IN)         :: HPROGRAM ! program calling surf. schemes
+INTEGER,          INTENT(IN)         :: KI       ! Surfex grid dimension
+END SUBROUTINE COUPLING_SURF_TOPD
+!     ######################
+!     ######################
+!     ######################
+SUBROUTINE PGD_TOPD(HPROGRAM)
+ CHARACTER(LEN=*),  INTENT(IN)     :: HPROGRAM    
+END SUBROUTINE PGD_TOPD
+!     ######################
+!     ######################
+!     ######################
+SUBROUTINE READ_NAMELISTS_TOPD(HPROGRAM)
+ CHARACTER(LEN=*),  INTENT(IN)     :: HPROGRAM    
+END SUBROUTINE READ_NAMELISTS_TOPD
+
 
diff --git a/src/SURFEX/budget_coupl_rout.F90 b/src/SURFEX/budget_coupl_rout.F90
deleted file mode 100644
index 12689fb02c9fd97912edaf41355b4e3751fb5466..0000000000000000000000000000000000000000
--- a/src/SURFEX/budget_coupl_rout.F90
+++ /dev/null
@@ -1,241 +0,0 @@
-
-!     ##########################
-      SUBROUTINE BUDGET_COUPL_ROUT(KNI,KFORC_STEP)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!    L. Bouilloud &  B. Vincendon	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original  03/2008 
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-USE MODD_BUDGET_COUPL_ROUT ! contains all useful variables XB_*
-USE MODD_TOPODYN,       ONLY : NNCAT, XQTOT, XTOPD_STEP, XQB_DR, XQB_RUN
-USE MODD_COUPLING_TOPD, ONLY : XRUNOFF_TOP, XRUN_TOROUT, XDR_TOROUT
-!
-USE MODD_SURF_PAR
-USE MODD_CSTS,            ONLY : XRHOLW
-!
-USE MODD_ISBA_n,          ONLY : XWG, XDG, XWR, CRUNOFF, XWGI, TSNOW
-USE MODD_DIAG_EVAP_ISBA_n,ONLY : XAVG_EVAPC, XAVG_LEGC, XAVG_RUNOFFC, XAVG_DRAINC,&
-                                 XAVG_DRAIN, XAVG_RUNOFF, XAVG_EVAP !ludo
-USE MODD_ISBA_GRID_n,     ONLY : XMESH_SIZE !ludo
-USE MODD_FORC_ATM,        ONLY : XSNOW       ,&! snow precipitation                    (kg/m2/s)
-                                 XRAIN         ! liquid precipitation                  (kg/m2/s)
-USE MODD_SURF_ATM_n,       ONLY:NR_NATURE
-!
-USE MODI_UNPACK_SAME_RANK
-USE MODI_PACK_SAME_RANK 
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-INTEGER, INTENT(IN)           :: KNI        ! expected physical size of full surface array
-INTEGER, INTENT(IN)           :: KFORC_STEP ! time step
-!
-!*      0.2    declarations of local variables
-REAL                          :: ZFACT0, ZFACT1, ZFACT2
-INTEGER                       :: JMESH,JCAT
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('BUDGET_COUPL_ROUT',0,ZHOOK_HANDLE)
-!
-!*       1.     Budget computation:
-!               ---------------
-!
-XB_RAIN(:)=XRAIN(:)  !kg/m2/s
-XB_SNOW(:)=XSNOW(:)  !kg/m2/s
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XWR(:,1), XB_WR)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_EVAPC, XB_EVAP)        
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_RUNOFFC,XB_RUNOFF_ISBA)
-IF (CRUNOFF=='TOPD') THEN
-  CALL UNPACK_SAME_RANK(NR_NATURE,XRUNOFF_TOP,XB_RUNOFF_TOPD)
-ELSE
-  XB_RUNOFF_TOPD = XB_RUNOFF_ISBA
-ENDIF
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINC, XB_DRAIN)
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,2,1)  ,XB_WG2)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,3,1)  ,XB_WG3)
-WHERE ( XB_WG2/=XUNDEF .AND. XB_DG2/=XUNDEF .AND. XB_WG3/=XUNDEF .AND. XB_DG3/=XUNDEF )
-  XB_WGTOT(:) = XB_WG2(:)*XB_DG2(:) + XB_WG3(:)*(XB_DG3(:)-XB_DG2(:)) !m3/m2
-ELSEWHERE
-  XB_WGTOT(:) = XUNDEF
-ENDWHERE
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XWGI(:,2,1) ,XB_WGI2)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWGI(:,3,1) ,XB_WGI3)
-WHERE ( XB_WGI2/=XUNDEF .AND. XB_DG2/=XUNDEF .AND. XB_WGI3/=XUNDEF .AND. XB_DG3/=XUNDEF )
-  XB_WGITOT(:) = XB_WGI2(:)*XB_DG2(:) + XB_WGI3(:)*(XB_DG3(:)-XB_DG2(:)) !m3/m2
-ELSEWHERE
-  XB_WGTOT(:) = XUNDEF
-ENDWHERE
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,TSNOW%WSNOW(:,1,1),XB_SWE1)
- CALL UNPACK_SAME_RANK(NR_NATURE,TSNOW%WSNOW(:,2,1),XB_SWE2)
- CALL UNPACK_SAME_RANK(NR_NATURE,TSNOW%WSNOW(:,3,1),XB_SWE3)
-XB_SWETOT(:) = XB_SWE1(:)+XB_SWE2(:)+XB_SWE3(:)
-!
-DO JCAT=1,NNCAT
-  XB_QTOT(JCAT) = SUM(XQTOT(JCAT,1:KFORC_STEP))
-  XB_QRUN(JCAT) = SUM(XQB_RUN(JCAT,1:KFORC_STEP))
-  XB_QDR(JCAT)  = SUM(XQB_DR(JCAT,1:KFORC_STEP))
-ENDDO
-!
-DO JCAT=1,NNCAT
-  !
-  DO JMESH=1,KNI
-    !
-    IF ( XB_DG2(JMESH)/=XUNDEF ) THEN
-!!    Water going in the system
-      ZFACT0 =  XB_ABV_BYMESH(JMESH,JCAT) * XB_MESH_SIZE(JMESH)
-      ZFACT1 =  ZFACT0 / XRHOLW
-      ZFACT2 =  XTOPD_STEP * ZFACT1
-      !
-!!      variable =1 : Rain      
-      XB_VAR_BV(KFORC_STEP,JCAT,1) = XB_VAR_BV(KFORC_STEP,JCAT,1) + XB_RAIN(JMESH) * ZFACT2
-!!      variable =2 : Snow
-      XB_VAR_BV(KFORC_STEP,JCAT,2) = XB_VAR_BV(KFORC_STEP,JCAT,2) + XB_SNOW(JMESH) * ZFACT2
-!!    Water going out of the system
-!!      variable =3 : Incerception 
-      XB_VAR_BV(KFORC_STEP,JCAT,3) = XB_VAR_BV(KFORC_STEP,JCAT,3) + (XB_WR(JMESH)-XB_WRM(JMESH)) * ZFACT1
-!!      variable =4 : Evaporation
-      XB_VAR_BV(KFORC_STEP,JCAT,4) = XB_VAR_BV(KFORC_STEP,JCAT,4) + (XB_EVAP(JMESH)-XB_EVAPM(JMESH)) * ZFACT1
-!!      variable =5 : Runoff
-      XB_VAR_BV(KFORC_STEP,JCAT,5) = XB_VAR_BV(KFORC_STEP,JCAT,5) + (XB_RUNOFF_TOPD(JMESH)-XB_RUNOFF_TOPDM(JMESH)) * ZFACT1
-!!      variable =6 : Drainage
-      XB_VAR_BV(KFORC_STEP,JCAT,6) = XB_VAR_BV(KFORC_STEP,JCAT,6) + (XB_DRAIN(JMESH)-XB_DRAINM(JMESH)) * ZFACT1
-!!      variable =7 : Variation of liquid water stocked in the ground
-      XB_VAR_BV(KFORC_STEP,JCAT,7) = XB_VAR_BV(KFORC_STEP,JCAT,7) + (XB_WGTOT(JMESH)-XB_WGTOTM(JMESH)) * ZFACT0   
-!!      variable =8 : Variation of solid water stocked in the ground
-      XB_VAR_BV(KFORC_STEP,JCAT,8) = XB_VAR_BV(KFORC_STEP,JCAT,8) + (XB_WGITOT(JMESH)-XB_WGITOTM(JMESH)) * ZFACT0 
-!!      variable =9 : Variation of melting snow
-      XB_VAR_BV(KFORC_STEP,JCAT,9) = XB_VAR_BV(KFORC_STEP,JCAT,9) + (XB_SWETOT(JMESH)-XB_SWETOTM(JMESH)) * ZFACT0 
-
- !! bilan hors BV en m3
-      ZFACT0 =  (1.-XB_ABV_BYMESH(JMESH,JCAT)) * XB_MESH_SIZE(JMESH)
-      ZFACT1 =  ZFACT0 / XRHOLW
-      ZFACT2 =  XTOPD_STEP * ZFACT1
-      !
-      XB_VAR_NOBV(KFORC_STEP,JCAT,1) = XB_VAR_NOBV(KFORC_STEP,JCAT,1) + XB_RAIN(JMESH) * ZFACT2
-      XB_VAR_NOBV(KFORC_STEP,JCAT,2) = XB_VAR_NOBV(KFORC_STEP,JCAT,2) + XB_SNOW(JMESH) * ZFACT2
-      XB_VAR_NOBV(KFORC_STEP,JCAT,3) = XB_VAR_NOBV(KFORC_STEP,JCAT,3) + (XB_WR(JMESH)-XB_WRM(JMESH)) * ZFACT1
-      XB_VAR_NOBV(KFORC_STEP,JCAT,4) = XB_VAR_NOBV(KFORC_STEP,JCAT,4) + (XB_EVAP(JMESH)-XB_EVAPM(JMESH)) * ZFACT1
-      XB_VAR_NOBV(KFORC_STEP,JCAT,5) = XB_VAR_NOBV(KFORC_STEP,JCAT,5) + (XB_RUNOFF_ISBA(JMESH)-XB_RUNOFF_ISBAM(JMESH)) * ZFACT1
-      XB_VAR_NOBV(KFORC_STEP,JCAT,6) = XB_VAR_NOBV(KFORC_STEP,JCAT,6) + (XB_DRAIN(JMESH)-XB_DRAINM(JMESH)) * ZFACT1
-      XB_VAR_NOBV(KFORC_STEP,JCAT,7) = XB_VAR_NOBV(KFORC_STEP,JCAT,7) + (XB_WGTOT(JMESH)-XB_WGTOTM(JMESH)) * ZFACT0   
-      XB_VAR_NOBV(KFORC_STEP,JCAT,8) = XB_VAR_NOBV(KFORC_STEP,JCAT,8) + (XB_WGITOT(JMESH)-XB_WGITOTM(JMESH)) * ZFACT0   
-      XB_VAR_NOBV(KFORC_STEP,JCAT,9) = XB_VAR_NOBV(KFORC_STEP,JCAT,9) + (XB_SWETOT(JMESH)-XB_SWETOTM(JMESH)) * ZFACT0 
-      !
-    ENDIF
-    !     
-  ENDDO
-  !
-  XB_VAR_BV(KFORC_STEP,JCAT,10) = SUM(XB_VAR_BV(KFORC_STEP,JCAT,1:2)) - SUM(XB_VAR_BV(KFORC_STEP,JCAT,3:9))
-  !
-  XB_VAR_NOBV(KFORC_STEP,JCAT,10) =  SUM(XB_VAR_NOBV(KFORC_STEP,JCAT,1:2)) - SUM(XB_VAR_NOBV(KFORC_STEP,JCAT,3:9))
-  !
-  XB_VAR_Q(KFORC_STEP,JCAT,1) = (XB_QTOT(JCAT)-XB_QTOTM(JCAT)) * XTOPD_STEP
-  XB_VAR_Q(KFORC_STEP,JCAT,2) = (XB_QRUN(JCAT)-XB_QRUNM(JCAT)) * XTOPD_STEP
-  XB_VAR_Q(KFORC_STEP,JCAT,3) = (XB_QDR(JCAT)- XB_QDRM(JCAT)) * XTOPD_STEP
-  XB_VAR_Q(KFORC_STEP,JCAT,4) = SUM(XRUN_TOROUT(JCAT,:)) * XTOPD_STEP
-  XB_VAR_Q(KFORC_STEP,JCAT,5) = SUM(XDR_TOROUT(JCAT,:)) * XTOPD_STEP
-  !
-ENDDO   
-!
-!bilan tot isba (m3)
-DO JMESH=1,KNI
-  !
-  IF (XB_DG2(JMESH)/=XUNDEF) THEN  
-    !
-    ZFACT0 =  XB_MESH_SIZE(JMESH)
-    ZFACT1 =  ZFACT0 / XRHOLW
-    ZFACT2 =  XTOPD_STEP * ZFACT1
-    !
-!!    Water going in the system
-!!      variable =1 : Rain
-    XB_VAR_TOT(KFORC_STEP,1) = XB_VAR_TOT(KFORC_STEP,1) + XB_RAIN(JMESH) * ZFACT2
-!!      variable =2 : Snow
-    XB_VAR_TOT(KFORC_STEP,2) = XB_VAR_TOT(KFORC_STEP,2) + XB_SNOW(JMESH) * ZFACT2
-!!    Water going out of the system
-!!      variable =3 : Incerception 
-    XB_VAR_TOT(KFORC_STEP,3) = XB_VAR_TOT(KFORC_STEP,3) + (XB_WR(JMESH)-XB_WRM(JMESH)) * ZFACT1
-!!      variable =4 : Evaporation
-    XB_VAR_TOT(KFORC_STEP,4) = XB_VAR_TOT(KFORC_STEP,4) + (XB_EVAP(JMESH)-XB_EVAPM(JMESH)) * ZFACT1
-!!      variable =5 : Runoff
-    XB_VAR_TOT(KFORC_STEP,5) = XB_VAR_TOT(KFORC_STEP,5) + &
-       (XB_RUNOFF_ISBA(JMESH)-XB_RUNOFF_ISBAM(JMESH)) * (1-XB_ATOP_BYMESH(JMESH)) * ZFACT1 + &
-       (XB_RUNOFF_TOPD(JMESH)-XB_RUNOFF_TOPDM(JMESH)) * XB_ATOP_BYMESH(JMESH) * ZFACT1
-!!      variable =6 : Drainage
-    XB_VAR_TOT(KFORC_STEP,6) = XB_VAR_TOT(KFORC_STEP,6) + (XB_DRAIN(JMESH)-XB_DRAINM(JMESH)) * ZFACT1
-!!      variable =7 : Variation of liquid water stocked in the ground
-    XB_VAR_TOT(KFORC_STEP,7) = XB_VAR_TOT(KFORC_STEP,7) + (XB_WGTOT(JMESH)-XB_WGTOTM(JMESH)) * ZFACT0  
-!!      variable =8 : Variation of solid water stocked in the ground
-    XB_VAR_TOT(KFORC_STEP,8) = XB_VAR_TOT(KFORC_STEP,8) + (XB_WGITOT(JMESH)-XB_WGITOTM(JMESH)) * ZFACT0  
-!!      variable =9 : Variation of melting snow
-    XB_VAR_TOT(KFORC_STEP,9) = XB_VAR_TOT(KFORC_STEP,9) + (XB_SWETOT(JMESH)-XB_SWETOTM(JMESH)) * ZFACT0  
-    !
-  ENDIF
-  !
-ENDDO !JMESH
-!
-XB_VAR_TOT(KFORC_STEP,10) = SUM(XB_VAR_TOT(KFORC_STEP,1:2)) - SUM(XB_VAR_TOT(KFORC_STEP,3:9)) 
-!
-XB_WRM  (:) = XB_WR(:)
-XB_EVAPM(:) = XB_EVAP(:)
-!
-XB_RUNOFF_TOPDM(:) = XB_RUNOFF_TOPD(:)
-XB_RUNOFF_ISBAM(:) = XB_RUNOFF_ISBA(:)
-XB_DRAINM      (:) = XB_DRAIN(:)
-!
-XB_WG2M   (:) = XB_WG2(:)
-XB_WG3M   (:) = XB_WG3(:)
-XB_WGTOTM (:) = XB_WGTOT(:)
-XB_WGI2M  (:) = XB_WGI2(:)
-XB_WGI3M  (:) = XB_WGI3(:)
-XB_WGITOTM(:) = XB_WGITOT(:)
-XB_SWE1M  (:) = XB_SWE1(:)
-XB_SWE2M  (:) = XB_SWE2(:)
-XB_SWE3M  (:) = XB_SWE3(:)
-XB_SWETOTM(:) = XB_SWETOT(:)
-!
-XB_QTOTM(:) = XB_QTOT(:)
-XB_QRUNM(:) = XB_QRUN(:) 
-XB_QDRM (:) = XB_QDR(:)
-!
-IF (LHOOK) CALL DR_HOOK('BUDGET_COUPL_ROUT',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE BUDGET_COUPL_ROUT
diff --git a/src/SURFEX/coupl_topd.F90 b/src/SURFEX/coupl_topd.F90
deleted file mode 100644
index 4671539d5779d3e4ed52472fa56daee51d2eebdf..0000000000000000000000000000000000000000
--- a/src/SURFEX/coupl_topd.F90
+++ /dev/null
@@ -1,383 +0,0 @@
-!-----------------------------------------------------------------
-!     #####################
-      SUBROUTINE COUPL_TOPD(HPROGRAM,HSTEP,KI,KSTEP)
-!     #####################
-!
-!!****  *COUPL_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!
-!     
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   15/10/2003
-!!      09/2007 : New organisation of exfiltration, computation of saturated
-!!                area, routing.
-!!                Soil ice content taken into account
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN,        ONLY : NNCAT, NMESHT, NNMC, XMPARA
-USE MODD_COUPLING_TOPD,  ONLY : XWG_FULL, XDTOPI, XKAC_PRE, XDTOPT, XWTOPT, XWSTOPT, XAS_NATURE,&
-                                  XKA_PRE, NMASKT, XWSUPSAT,&
-                                  XRUNOFF_TOP, XATOP, XWFCTOPI, NNPIX,&
-                                  XFRAC_D2, XFRAC_D3, XWSTOPI, XDMAXFC, XWFCTOPT, XWGI_FULL,&
-                                  NFREQ_MAPS_ASAT, XAVG_RUNOFFCM, XAVG_DRAINCM
-!
-USE MODD_ISBA_n,           ONLY : XWG, XDG, XWGI
-USE MODD_CSTS,             ONLY : XRHOLW, XRHOLI
-USE MODD_DIAG_EVAP_ISBA_n, ONLY : XAVG_RUNOFFC, XAVG_DRAINC
-USE MODD_SURF_ATM_n,       ONLY : NR_NATURE
-USE MODD_SURF_ATM_GRID_n,  ONLY : XMESH_SIZE
-USE MODD_SURF_PAR,         ONLY : XUNDEF, NUNDEF
-USE MODD_ISBA_PAR,         ONLY : XWGMIN
-!
-USE MODI_GET_LUOUT
-USE MODI_UNPACK_SAME_RANK
-USE MODI_PACK_SAME_RANK
-USE MODI_ISBA_TO_TOPD
-USE MODI_RECHARGE_SURF_TOPD
-USE MODI_TOPODYN_LAT
-USE MODI_SAT_AREA_FRAC
-USE MODI_TOPD_TO_ISBA
-USE MODI_DIAG_ISBA_TO_ROUT
-USE MODI_ISBA_TO_TOPDSAT
-USE MODI_ROUTING
-USE MODI_OPEN_FILE
-USE MODI_WRITE_FILE_ISBAMAP
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
- CHARACTER(LEN=*), INTENT(IN) :: HSTEP  ! atmospheric loop index
-INTEGER, INTENT(IN)          :: KI    ! Grid dimensions
-INTEGER, INTENT(IN)          :: KSTEP ! current time step 
-!
-!*      0.2    declarations of local variables
-!
-REAL, DIMENSION(SIZE(XRUNOFF_TOP,1)) :: ZWG2_TMP,ZWG3_TMP
-REAL, DIMENSION(SIZE(XRUNOFF_TOP,1)) :: ZWSAT_NAT
-REAL, DIMENSION(NNCAT,NMESHT) :: ZRT             ! recharge on TOP-LAT grid (m)
-REAL, DIMENSION(NNCAT,NMESHT) :: ZDEFT           ! local deficits on TOPODYN grid (m)
-REAL, DIMENSION(NNCAT,NMESHT) :: ZRI_WGIT        ! water changing of phase on TOPMODEL grid
-REAL, DIMENSION(NNCAT,NMESHT) :: ZRUNOFF_TOPD    ! Runoff on the Topodyn grid (m3/s)
-REAL, DIMENSION(NNCAT,NMESHT) :: ZDRAIN_TOPD     ! Drainage from Isba on Topodyn grid (m3/s)
-REAL, DIMENSION(NNCAT,NMESHT) :: ZKAPPA          ! topographic index
-REAL, DIMENSION(NNCAT)        :: ZKAPPAC         ! critical topographic index
-REAL, DIMENSION(KI)           :: ZRI             ! recharge on ISBA grid (m)
-REAL, DIMENSION(KI)           :: ZRI_WGI         ! water changing of phase on ISBA grid
-REAL, DIMENSION(KI)           :: Z_WSTOPI, Z_WFCTOPI
-REAL, DIMENSION(KI)           :: ZWM             ! Water content on SurfEx grid after the previous topodyn time step
-REAL, DIMENSION(KI)           :: ZWGIM           ! soil ice content at previous time step
-REAL, DIMENSION(KI)           :: ZRUNOFFC_FULL   ! Cumulated runoff from isba on the full domain (kg/m2)
-REAL, DIMENSION(KI)           :: ZRUNOFFC_FULLM  ! Cumulated runoff from isba on the full domain (kg/m2) at t-dt
-REAL, DIMENSION(KI)           :: ZRUNOFF_ISBA    ! Runoff from Isba (kg/m2)
-REAL, DIMENSION(KI)           :: ZDRAINC_FULL    ! Cumulated drainage from Isba on the full domain (kg/m2)
-REAL, DIMENSION(KI)           :: ZDRAINC_FULLM   ! Cumulated drainage from Isba on the full domain (kg/m2) at t-dt
-REAL, DIMENSION(KI)           :: ZDRAIN_ISBA     ! Drainage from Isba (m3/s)
-REAL, DIMENSION(KI)           :: ZDG_FULL
-REAL, DIMENSION(KI)           :: ZWG2_FULL, ZWG3_FULL, ZDG2_FULL, ZDG3_FULL
-REAL, DIMENSION(KI)           :: ZWGI_FULL
-REAL, DIMENSION(KI)           :: ZAS             ! Saturated area fraction for each Isba meshes
-REAL, DIMENSION(NNCAT)        :: Z_DW1,Z_DW2     ! Wsat-Wfc to actualise M in fonction of WI
-REAL                          :: ZAVG_MESH_SIZE
-LOGICAL, DIMENSION(NNCAT)     :: GTOPD           ! logical variable = true if topodyn_lat runs
-INTEGER                       :: JJ,JI           ! loop control 
-INTEGER                       :: IUNIT           ! unit number of results files
-INTEGER                       :: ILUOUT          ! unit number of listing file
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('COUPL_TOPD',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-!*       0.     Initialization:
-!               ---------------
-!
-ZWM(:) = XWG_FULL(:)
-ZWGIM(:) = XWGI_FULL(:)
-ZWG2_TMP(:) = XWG(:,2,1)
-ZWG3_TMP(:) = XWG(:,3,1)
-!
-!*       1.     ISBA => TOPODYN
-!               ---------------
-!*       1.1    Computation of the useful depth and water for lateral transfers
-!               -----------------------------------
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,2,1),ZDG2_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,3,1),ZDG3_FULL)
-WHERE ( ZDG2_FULL/=XUNDEF )
-  ZDG_FULL = XFRAC_D2*ZDG2_FULL + XFRAC_D3*(ZDG3_FULL-ZDG2_FULL)
-ELSEWHERE
-  ZDG_FULL = XUNDEF
-END WHERE
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,2,1),ZWG2_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,3,1),ZWG3_FULL)
-WHERE ( ZDG_FULL/=XUNDEF .AND. ZDG_FULL/=0. )
-  XWG_FULL = XFRAC_D2*(ZDG2_FULL/ZDG_FULL)*ZWG2_FULL + XFRAC_D3*((ZDG3_FULL-ZDG2_FULL)/ZDG_FULL)*ZWG3_FULL
-ELSEWHERE
-  XWG_FULL = XUNDEF
-END WHERE
-!
-!ludo prise en compte glace (pas de glace dans 3e couche)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWGI(:,2,1),ZWGI_FULL)
-WHERE ( ZWGI_FULL/=XUNDEF .AND. XFRAC_D2>0 .AND. ZDG_FULL/=0. )
-  XWGI_FULL = XFRAC_D2*(ZDG2_FULL/ZDG_FULL)*ZWGI_FULL
-ELSEWHERE
-  XWGI_FULL = XUNDEF
-END WHERE
-!
-!
-WHERE ( XDTOPI/=XUNDEF ) 
-  ZRI_WGI = ( XWGI_FULL - ZWGIM ) !m3/m3
-ELSEWHERE
-  ZRI_WGI = 0.0
-END WHERE
-!
- CALL ISBA_TO_TOPD(ZRI_WGI,ZRI_WGIT)
-!
-!*       1.2    Water recharge 
-!               ---------------
-! Topodyn uses :
-! - a water recharge = water added since last time step to compute hydrological similarity indexes
-! - the total water content to compute a deficit
-!
-! This recharge is computed without regarding the changing of phase of water
-! and the lateral transfers are performed regarding wsat et Wfc of last time step
-!
-WHERE ( XDTOPI/=XUNDEF )
-  ZRI = ( (XWG_FULL - ZWM) + ZRI_WGI ) * XDTOPI
-ELSEWHERE
-  ZRI = 0.0
-ENDWHERE
-!
-! The water recharge on ISBA grid is computed on TOPMODEL grid
- CALL RECHARGE_SURF_TOPD(ZRI,ZRT,KI)
-!
-!*       2.     Lateral distribution
-!               --------------------
-!*       2.1    Computation of local deficits on TOPODYN grid
-!               ----------------------------------------
-!
- CALL TOPODYN_LAT(ZRT(:,:),ZDEFT(:,:),ZKAPPA(:,:),ZKAPPAC(:),GTOPD)
-!
-!*       2.2    Computation of contributive area on ISBA grid
-!               ----------------------------------------
-!
- CALL SAT_AREA_FRAC(ZDEFT,ZAS)
-!
- CALL PACK_SAME_RANK(NR_NATURE,ZAS,XAS_NATURE)
-!
-!*       3.    Deficit (m) -> water storage (m3/m3) and changing of phase
-!               ------------------------------------
-!ancien Wsat-Wfc pour actualisation param M
-Z_DW1(:)=0
-DO JJ=1,NNCAT
-  Z_DW1(JJ) = SUM( XWSTOPT(JJ,:)-XWFCTOPT(JJ,:), MASK=XWSTOPT(JJ,:)/=XUNDEF ) / NNMC(JJ)
-ENDDO
-!
-DO JJ=1,NNCAT
-  WHERE ( XDTOPT(JJ,:)/=XUNDEF .AND. XDTOPT(JJ,:)/=0. )
-    XWTOPT(JJ,:) = XWSTOPT(JJ,:) - ( ZDEFT(JJ,:) / XDTOPT(JJ,:) )      
-    !changing phase
-    XWTOPT(JJ,:) = XWTOPT(JJ,:) - ZRI_WGIT(JJ,:)
-  END WHERE
-ENDDO
-WHERE (XWTOPT > XWSTOPT ) XWTOPT = XWSTOPT
-!
-!actualisation Wsat, Wfc et Dmax pour prochain pas de temps
-WHERE ( ZWG2_FULL/=XUNDEF .AND. XWSTOPI/=0. )
-  Z_WSTOPI  = XWSTOPI - XWGI_FULL
-  Z_WFCTOPI = XWFCTOPI * Z_WSTOPI / XWSTOPI
-END WHERE
- CALL ISBA_TO_TOPD(Z_WSTOPI,XWSTOPT)
- CALL ISBA_TO_TOPD(Z_WFCTOPI,XWFCTOPT)
-!
-!ludo test empeche erreur num chgt phase
-WHERE ( ABS(XWSTOPT-XWTOPT) < 0.0000000001 ) XWSTOPT = XWTOPT
-!
-WHERE ( XWTOPT>XWSTOPT ) XWTOPT = XWSTOPT
-!
-WHERE ( XWFCTOPT/= XUNDEF ) XDMAXFC = (XWSTOPT - XWFCTOPT) * XDTOPT ! (m)
-!
-!actualisation M-> M=M*rapport Wsat-Wfc
-Z_DW2(:)=0
-DO JJ=1,NNCAT
-  Z_DW2(JJ) = SUM( XWSTOPT(JJ,:)-XWFCTOPT(JJ,:), MASK=XWSTOPT(JJ,:)/=XUNDEF) / NNMC(JJ)
-ENDDO
-!
-XMPARA(:) = XMPARA(:) * Z_DW2(:)/Z_DW1(:)
-!
-!*       4.     TOPODYN => ISBA
-!               ---------------
-!*       4.1    Calculation of water storage on ISBA grid
-!               -----------------------------------------
-!
- CALL TOPD_TO_ISBA(KI,KSTEP,GTOPD)
- CALL PACK_SAME_RANK(NR_NATURE, (1-XFRAC_D2)*ZWG2_FULL + XFRAC_D2*XWG_FULL, XWG(:,2,1))
- CALL PACK_SAME_RANK(NR_NATURE, (1-XFRAC_D3)*ZWG3_FULL + XFRAC_D3*XWG_FULL, XWG(:,3,1))
-!
-!*       4.2    Budget correction
-!  -----------------------------------------
-!
-ZAVG_MESH_SIZE = SUM(XMESH_SIZE(:)) / SIZE(XMESH_SIZE,1)
-!
- CALL BUDG_WG(ZWG2_TMP(:),XWG(:,2,1),XDG(:,2,1),XMESH_SIZE,ZAVG_MESH_SIZE)
-!   
- CALL BUDG_WG(ZWG3_TMP(:),XWG(:,3,1),XDG(:,3,1)-XDG(:,2,1),XMESH_SIZE,ZAVG_MESH_SIZE)
-!
-!*      5.0    Total discharge
-!              ---------------
-!
-!*      5.1    Total water for runoff on TOPODYN grid
-!              ---------------------------------------
-!
- CALL PACK_SAME_RANK(NR_NATURE,XWSUPSAT,ZWSAT_NAT)
-WHERE ( ZWSAT_NAT(:)<10E-10 ) 
-  ZWSAT_NAT(:) = 0.
-ELSEWHERE( ZWSAT_NAT(:)/=XUNDEF ) 
-  XRUNOFF_TOP(:) = XRUNOFF_TOP(:) + ZWSAT_NAT(:)*XRHOLW*XDG(:,2,1)
-ENDWHERE
-!
-!write(*,*) ' ds coupl_topd ZWSAT_NAT',SUM(ZWSAT_NAT,MASK=ZWSAT_NAT(:)/=XUNDEF)!
- CALL UNPACK_SAME_RANK(NR_NATURE,XRUNOFF_TOP,ZRUNOFFC_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_RUNOFFCM,ZRUNOFFC_FULLM)
-!
- CALL DIAG_ISBA_TO_ROUT(ZRUNOFFC_FULL,ZRUNOFFC_FULLM,ZRUNOFF_ISBA)
-!
-XAVG_RUNOFFCM(:) = XRUNOFF_TOP(:)
-!
-WHERE (ZRUNOFF_ISBA==XUNDEF) ZRUNOFF_ISBA = 0.
-!
-ZRUNOFF_TOPD(:,:) = 0
-!
- CALL ISBA_TO_TOPDSAT(XKA_PRE,XKAC_PRE,KI,ZRUNOFF_ISBA,ZRUNOFF_TOPD)
-!
-!
-!*      5.2    Total water for drainage on TOPODYN grid
-!              ----------------------------------------
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINC*XATOP,ZDRAINC_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINCM*XATOP,ZDRAINC_FULLM)
-!
- CALL DIAG_ISBA_TO_ROUT(ZDRAINC_FULL,ZDRAINC_FULLM,ZDRAIN_ISBA)
-!
-WHERE (ZDRAIN_ISBA==XUNDEF) ZDRAIN_ISBA=0.
-!
-XAVG_DRAINCM(:)  = XAVG_DRAINC(:)
-!
-ZDRAIN_TOPD(:,:) = 0.0
-!
- CALL ISBA_TO_TOPD(ZDRAIN_ISBA,ZDRAIN_TOPD)
-!
-DO JJ=1,NNCAT
-  DO JI=1,NNMC(JJ)
-    IF (NMASKT(JJ,JI)/=NUNDEF) &
-      ZDRAIN_TOPD(JJ,JI) = ZDRAIN_TOPD(JJ,JI) / NNPIX(NMASKT(JJ,JI))
-  ENDDO
-ENDDO   
-!
-!*      6    Routing (runoff + drainage + exfiltration)
-!
- CALL ROUTING(ZRUNOFF_TOPD,ZDRAIN_TOPD,KSTEP)
-!
-XKA_PRE(:,:) = ZKAPPA(:,:)
-XKAC_PRE(:) = ZKAPPAC(:)
-!
-!*      7.0    Writing results in map files
-!              ----------------------------
-!
-IF (NFREQ_MAPS_ASAT/=0.AND.MOD(KSTEP,NFREQ_MAPS_ASAT)==0) THEN
-  CALL OPEN_FILE('ASCII ',IUNIT,HFILE='carte_surfcont'//HSTEP,HFORM='FORMATTED',HACTION='WRITE')
-  CALL WRITE_FILE_ISBAMAP(IUNIT,ZAS,KI)
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-ENDIF
-!
-IF (LHOOK) CALL DR_HOOK('COUPL_TOPD',1,ZHOOK_HANDLE)
-!
-CONTAINS
-!
-SUBROUTINE BUDG_WG(PWGM,PWG,PDG,PMESH_SIZE,PAVG_MESH_SIZE)
-!
-REAL, DIMENSION(:), INTENT(IN) :: PWGM
-REAL, DIMENSION(:), INTENT(INOUT) :: PWG
-REAL, DIMENSION(:), INTENT(IN) :: PDG
-REAL, DIMENSION(:), INTENT(IN) :: PMESH_SIZE
-REAL, INTENT(IN) :: PAVG_MESH_SIZE
-!
-REAL :: ZSTOCK_WGM, ZSTOCK_WG
-REAL :: ZAVG_DG, ZBUDG_WG
-REAL :: ZTMP, ZTMP2
-INTEGER :: JMESH
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-IF (LHOOK) CALL DR_HOOK('COUPL_TOPD:BUDG_WG',0,ZHOOK_HANDLE)
-!
-ZSTOCK_WGM = SUM( PWGM(:)*PDG(:)*PMESH_SIZE(:) )    ! water stocked in the ground (m3)
-ZSTOCK_WG  = SUM( PWG (:)*PDG(:)*PMESH_SIZE(:) )    ! water stocked in the layer2 (m3)
-ZAVG_DG = SUM(PDG(:)) / SIZE(PDG)
-ZBUDG_WG = ( ZSTOCK_WG - ZSTOCK_WGM )/ ZAVG_DG / PAVG_MESH_SIZE
-!ZTMP2=0.
-!DO JMESH=1,SIZE(PWG)             
-! IF ((PWG(JMESH)/=XUNDEF).AND.(PWGM(JMESH)/=XUNDEF)) THEN
-!   IF(PWG(JMESH)/=PWGM(JMESH)) ZTMP2=ZTMP2+1.
-! ENDIF
-!ENDDO
-!
-ZTMP  = COUNT( PWG(:)/=PWGM(:) )
-ZTMP2 = COUNT( PWG(:)/=PWGM(:) .AND. PWG(:)>XWGMIN+(ZBUDG_WG/ZTMP) )
-!   
-DO JMESH=1,SIZE(PWG)
-  IF( PWG(JMESH)/=PWGM(JMESH) .AND. ZTMP2/=0. .AND. PWG(JMESH)>XWGMIN+(ZBUDG_WG/ZTMP) ) THEN
-    !IF( PWG(JMESH)/=PWGM(JMESH) .AND. (ZTMP2/=0. ) THEN
-    PWG(JMESH) = PWG(JMESH) - (ZBUDG_WG/ZTMP2)
-    !PWG(JMESH) = MAX(PWG(JMESH) - (ZBUDG_WG/ZTMP2),XWGMIN)
-  ENDIF
-ENDDO
-!write(*,*) KSTEP,COUNT(PWG(:)<XWGMIN),COUNT(PWG(:)<0.)
-!
-IF (LHOOK) CALL DR_HOOK('COUPL_TOPD:BUDG_WG',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE BUDG_WG
-!
-END SUBROUTINE COUPL_TOPD
diff --git a/src/SURFEX/coupling_surf_topd.F90 b/src/SURFEX/coupling_surf_topd.F90
deleted file mode 100644
index b1d82b3ed6835a1b335685319f581bbc41d9c0a4..0000000000000000000000000000000000000000
--- a/src/SURFEX/coupling_surf_topd.F90
+++ /dev/null
@@ -1,98 +0,0 @@
-!###################################################################
-SUBROUTINE COUPLING_SURF_TOPD (HPROGRAM,KI)
-!###################################################################
-!
-!!****  *COUPLING_SURF_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!!   
-!!    Driver for the coupling between SURFEX and TOPODYN
-!!      
-!!    REFERENCE
-!!    ---------
-!!    *COUPLING_SURF_TRIP from B. Decharme
-!!      
-!!    AUTHOR
-!!    ------
-!!	B. Vincendon    
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    07/06/11 
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODI_GET_LUOUT
-USE MODI_COUPL_TOPD
-USE MODI_ROUT_DATA_ISBA
-USE MODI_BUDGET_COUPL_ROUT
-USE MODI_WRITE_DISCHARGE_FILE
-USE MODI_WRITE_BUDGET_COUPL_ROUT
-USE MODI_PREP_RESTART_COUPL_TOPD
-!
-USE MODD_TOPODYN,       ONLY : XQTOT, NNB_TOPD_STEP, XQB_RUN, XQB_DR
-USE MODD_COUPLING_TOPD, ONLY : LCOUPL_TOPD, LBUDGET_TOPD, NNB_TOPD, LTOPD_STEP, NTOPD_STEP, &
-                                 NYEAR,NMONTH,NDAY,NH,NM
-!
-USE MODD_ISBA_n,          ONLY : CRUNOFF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=6), INTENT(IN)         :: HPROGRAM ! program calling surf. schemes
-INTEGER,          INTENT(IN)         :: KI       ! Surfex grid dimension
-                                                 ! in a forcing iteration
-!
-!*      0.2    declarations of local variables
-!
- CHARACTER(LEN=3)              :: YSTEP    ! time stepsurf_tmp/off
-INTEGER                       :: ILUOUT   ! unit of output listing file
-INTEGER                       :: JJ       ! loop control
-!
-REAL, DIMENSION(KI)           :: ZDG_FULL
-REAL, DIMENSION(KI)           :: ZWG2_FULL,ZWG3_FULL,ZDG2_FULL,ZDG3_FULL
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('COUPLING_SURF_TOPD',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-IF ( .NOT.LCOUPL_TOPD ) THEN
-  IF (LHOOK) CALL DR_HOOK('COUPLING_SURF_TOPD',1,ZHOOK_HANDLE)
-  RETURN
-ENDIF
-  !
-IF ( LTOPD_STEP ) THEN
-  !
-  ! * 1. Calling coupling or routing
-  !
-  IF (NTOPD_STEP<10) THEN
-    WRITE(YSTEP,'(I1)') NTOPD_STEP
-  ELSEIF (NTOPD_STEP < 100) THEN
-    WRITE(YSTEP,'(I2)') NTOPD_STEP
-  ELSE
-    WRITE(YSTEP,'(I3)') NTOPD_STEP
-  ENDIF
-  !
-  write(*,*) 'pas de temps coupl ',YSTEP
-  !
-  IF (CRUNOFF=='TOPD') THEN
-    CALL COUPL_TOPD(HPROGRAM,YSTEP,KI,NTOPD_STEP)
-  ELSE
-    CALL ROUT_DATA_ISBA(HPROGRAM,KI,NTOPD_STEP)
-  ENDIF
-  !
-  IF (LBUDGET_TOPD) CALL BUDGET_COUPL_ROUT(KI,NTOPD_STEP)
-  !
-ENDIF! (LCOUPL_TOPD.AND......
-!
-IF (LHOOK) CALL DR_HOOK('COUPLING_SURF_TOPD',1,ZHOOK_HANDLE)
-!-------------------------------------------------------------------------------
-END SUBROUTINE COUPLING_SURF_TOPD
diff --git a/src/SURFEX/diag_isba_to_rout.F90 b/src/SURFEX/diag_isba_to_rout.F90
deleted file mode 100644
index d702eef0bc2c21e08f73469930d09b319538fbaf..0000000000000000000000000000000000000000
--- a/src/SURFEX/diag_isba_to_rout.F90
+++ /dev/null
@@ -1,93 +0,0 @@
-!-----------------------------------------------------------------
-!     ############################
-      SUBROUTINE DIAG_ISBA_TO_ROUT(PVARC,PVARCP,PVARROUT)
-!     ############################
-!
-!!****  *DIAG_ISBA_TO_ROUT*  
-!!
-!!    PURPOSE
-!!    -------
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   10/11/2006
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_SURF_ATM_GRID_n, ONLY: XMESH_SIZE
-USE MODD_SURF_PAR,        ONLY: XUNDEF
-USE MODD_CSTS,            ONLY: XRHOLW
-USE MODD_TOPODYN, ONLY : XTOPD_STEP
-!
-USE MODI_ABOR1_SFX
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-!
-REAL,DIMENSION(:),INTENT(IN)        :: PVARC       ! Current time step cumulated diagnostic from SurfEx
-REAL,DIMENSION(:),INTENT(IN)        :: PVARCP      ! Previous time step cumulated diagnostic from SurfEx 
-REAL,DIMENSION(:),INTENT(OUT)       :: PVARROUT    ! Not cumulated diagnostic (m3/s)
-!
-!*      0.2    declarations of local variables
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('DIAG_ISBA_TO_ROUT',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-PVARROUT=XUNDEF
-!
-IF ( SIZE(PVARC,1)==SIZE(PVARCP,1) ) THEN
-  !
-  WHERE ( PVARC/=XUNDEF )
-    PVARROUT = PVARC - PVARCP
-    PVARROUT = PVARROUT / XTOPD_STEP
-    PVARROUT = PVARROUT * XMESH_SIZE / XRHOLW
-  ENDWHERE
-  !
-ELSE 
-  !
-  WRITE(*,*) 'Pb with diagnostic to rout'
-  CALL ABOR1_SFX("DIAG_ISBA_TO_ROUT: PB WITH DIAGNOSTIC TO ROUT ")
-  !
-ENDIF
-!
-WHERE (PVARROUT<0.) PVARROUT = 0.
-!
-IF (LHOOK) CALL DR_HOOK('DIAG_ISBA_TO_ROUT',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE DIAG_ISBA_TO_ROUT
diff --git a/src/SURFEX/flowdown.F90 b/src/SURFEX/flowdown.F90
deleted file mode 100644
index 71c6533ffb5010bf288f1445a57850fdc7f3f4f2..0000000000000000000000000000000000000000
--- a/src/SURFEX/flowdown.F90
+++ /dev/null
@@ -1,84 +0,0 @@
-!-----------------------------------------------------------------
-!     ###################
-      SUBROUTINE FLOWDOWN(KNMC,PVAR,PCONN,KLINE)
-!     ###################
-!
-!!****  *FLOWDOWN*
-!
-!!    PURPOSE
-!!    -------
-! to propagate data between pixels of a catchment in function of its topography
-!
-!
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* CNRM / Meteo-France *
-!!      G-M Saulnier    * LTHE *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   14/01/2005
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1     declarations of arguments
-!
-INTEGER, INTENT(IN) :: KNMC  ! catchment grid points number
-REAL, DIMENSION(:), INTENT(INOUT)   :: PVAR  ! variable to propagate
-REAL, DIMENSION(:,:), INTENT(IN)    :: PCONN ! catchment grid points connections
-INTEGER, DIMENSION(:), INTENT(IN)   :: KLINE ! 
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                  :: JJ, JI ! work variables
-INTEGER                  :: JNUP  ! number of upslope pixels
-INTEGER                  :: JCOL  ! third index of the pixel in the array XCONN
-INTEGER                  :: JREF  ! index of the upslope pixel in the topo domain
-REAL                     :: ZFAC  ! propagation factor between this pixel and the
-                                  ! upslope one
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('FLOWDOWN',0,ZHOOK_HANDLE)
-!
-DO JJ=1,KNMC
-   JNUP = INT(PCONN(JJ,4))
-   DO JI=1,JNUP
-      JCOL = ((JI-1)*2) + 5
-      JREF = INT(PCONN(JJ,JCOL))
-      ZFAC = PCONN(JJ,JCOL+1)
-      PVAR(JJ) = PVAR(JJ) + PVAR(KLINE(JREF)) * ZFAC
-   ENDDO
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('FLOWDOWN',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE FLOWDOWN
diff --git a/src/SURFEX/init_budget_coupl_rout.F90 b/src/SURFEX/init_budget_coupl_rout.F90
deleted file mode 100644
index 823a4162f947f425f2332222266e90df80b4d96e..0000000000000000000000000000000000000000
--- a/src/SURFEX/init_budget_coupl_rout.F90
+++ /dev/null
@@ -1,234 +0,0 @@
-!     ##########################
-      SUBROUTINE INIT_BUDGET_COUPL_ROUT(KNI)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!    Initialise varriables usefull for budget computation    
-!     
-!!**  METHOD
-!!    ------
-!!    Terms of the budget on all the domain
-!!    XB_VAR_TOT(forcing time step,variable)
-!!    Water going in the system
-!!      variable =1 : Rain
-!!      variable =2 : Snow
-!!    Water going out of the system
-!!      variable =3 : Incerception 
-!!      variable =4 : Evaporation
-!!      variable =5 : Runoff
-!!      variable =6 : Drainage
-!!      variable =7 : Variation of liquid water stocked in the ground
-!!      variable =8 : Variation of solid water stocked in the ground
-!!      variable =9 : Variation of melting snow
-!!    Budget 
-!!      variable =10:  Water going in the system- Water going out of the system
-!!    
-!!    Terms of the budget on a given catchment
-!!    XB_VAR_BV(forcing time step,catchment,variable)
-!!    XB_VAR_NOBV(forcing time step,catchment,variable)
-!!    
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!    L. Bouilloud &  B. Vincendon	* Meteo-France *
-
-!!    
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original  03/2008 
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-      ! declarative modules
-USE MODD_BUDGET_COUPL_ROUT ! contains all useful variables XB_*
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, NNMC, XQTOT,&
-                                 NNB_TOPD_STEP, XDXT, XQB_DR, XQB_RUN
-USE MODD_COUPLING_TOPD, ONLY : XAS_NATURE,  NNB_TOPD,&
-                                 XRUNOFF_TOP, XATOP, XBV_IN_MESH
-USE MODD_ISBA_n,          ONLY : XWG, XDG,  XWR, CRUNOFF, XWGI, TSNOW
-USE MODD_DIAG_EVAP_ISBA_n,ONLY : XAVG_EVAPC, XAVG_LEGC, XAVG_RUNOFFC ,XAVG_DRAINC,&
-                                 XAVG_DRAIN, XAVG_RUNOFF, XAVG_EVAP !ludo
-USE MODD_ISBA_GRID_n,     ONLY : XMESH_SIZE !ludo
-USE MODD_SURF_ATM_n,      ONLY : NR_NATURE
-!
-USE MODD_SURF_PAR,         ONLY:XUNDEF
-!
-USE MODI_UNPACK_SAME_RANK
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN)           :: KNI      ! expected physical size of full surface array
-!
-!*      0.2    declarations of local variables
-!
-INTEGER    :: JJ,JWRK2
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('INIT_BUDGET_COUPL_ROUT',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
-!write(*,*) 'In init_budget_coupl_rout KNI=',KNI
-ALLOCATE(YB_VAR(10)) 
-YB_VAR(1)='RAIN  '
-YB_VAR(2)='SNOW  '
-YB_VAR(3)='INTERC'
-YB_VAR(4)='EVATRA'
-YB_VAR(5)='RUNOFF'
-YB_VAR(6)='DRAINA'
-YB_VAR(7)='DSTOWG'
-YB_VAR(8)='DSTOWI'
-YB_VAR(9)='DSTOSW'
-YB_VAR(10)='BUDGET'
-!   
-!
-ALLOCATE(XB_RAIN(KNI))
-ALLOCATE(XB_SNOW(KNI))
-!
-ALLOCATE(XB_WR(KNI))
-ALLOCATE(XB_EVAP(KNI))
-ALLOCATE(XB_RUNOFF_ISBA(KNI))
-ALLOCATE(XB_DRAIN(KNI))
-ALLOCATE(XB_WG2(KNI))
-ALLOCATE(XB_WG3(KNI))
-ALLOCATE(XB_WGTOT(KNI))
-ALLOCATE(XB_WGI2(KNI))
-ALLOCATE(XB_WGI3(KNI))
-ALLOCATE(XB_WGITOT(KNI))
-ALLOCATE(XB_SWE1(KNI))
-ALLOCATE(XB_SWE2(KNI))
-ALLOCATE(XB_SWE3(KNI))
-ALLOCATE(XB_SWETOT(KNI))
-!
-ALLOCATE(XB_WRM(KNI))
-ALLOCATE(XB_EVAPM(KNI))
-ALLOCATE(XB_DRAINM(KNI))
-ALLOCATE(XB_RUNOFF_ISBAM(KNI))
-ALLOCATE(XB_WG2M(KNI))
-ALLOCATE(XB_WG3M(KNI))
-ALLOCATE(XB_WGTOTM(KNI))
-ALLOCATE(XB_WGI2M(KNI))
-ALLOCATE(XB_WGI3M(KNI))
-ALLOCATE(XB_WGITOTM(KNI))
-ALLOCATE(XB_SWE1M(KNI))
-ALLOCATE(XB_SWE2M(KNI))
-ALLOCATE(XB_SWE3M(KNI))   
-ALLOCATE(XB_SWETOTM(KNI))
-!
-ALLOCATE(XB_MESH_SIZE(KNI))
-ALLOCATE(XB_DG2(KNI))
-ALLOCATE(XB_DG3(KNI))
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XWR(:,1),XB_WRM)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_EVAPC,XB_EVAPM)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_RUNOFFC,XB_RUNOFF_ISBAM)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINC,XB_DRAINM)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,2,1),XB_WG2M)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,3,1),XB_WG3M)
- CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,2,1),XB_DG2)
- CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,3,1),XB_DG3)
-!
-WHERE ( XB_WG2M/=XUNDEF .AND. XB_DG2/=XUNDEF .AND. XB_WG3M/=XUNDEF .AND. XB_DG3/=XUNDEF )
-  XB_WGTOTM(:) = XB_WG2M(:)*XB_DG2(:) + XB_WG3M(:)*(XB_DG3(:)-XB_DG2(:)) !m3/m2
-ELSEWHERE
-  XB_WGTOTM(:) = XUNDEF
-ENDWHERE
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XWGI(:,2,1),XB_WGI2M)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWGI(:,3,1),XB_WGI3M)
-WHERE ((XB_WGI2M/=XUNDEF).AND.(XB_DG2/=XUNDEF).AND.(XB_WGI3M/=XUNDEF).AND.(XB_DG3/=XUNDEF))
-  XB_WGITOTM(:) = XB_WGI2M(:)*XB_DG2(:) + XB_WGI3M(:)*(XB_DG3(:)-XB_DG2(:)) !m3/m2
-ELSEWHERE
-  XB_WGITOTM(:) = XUNDEF
-ENDWHERE
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,TSNOW%WSNOW(:,1,1),XB_SWE1M)
- CALL UNPACK_SAME_RANK(NR_NATURE,TSNOW%WSNOW(:,2,1),XB_SWE2M)
- CALL UNPACK_SAME_RANK(NR_NATURE,TSNOW%WSNOW(:,3,1),XB_SWE3M)
-XB_SWETOTM(:) = XB_SWE1M(:)+XB_SWE2M(:)+XB_SWE3M(:)
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XMESH_SIZE,XB_MESH_SIZE)
-!
-ALLOCATE(XB_ABV_BYMESH(KNI,NNCAT))
-DO JJ=1,KNI
-  XB_ABV_BYMESH(JJ,:) = XBV_IN_MESH(JJ,:)/XB_MESH_SIZE(JJ) !*NNMC(:)*XDXT(:)**2 ! 
-  XB_ABV_BYMESH(JJ,:) = MIN(1.,XB_ABV_BYMESH(JJ,:))      
-ENDDO
-! 
-ALLOCATE(XB_VAR_BV(NNB_TOPD_STEP,NNCAT,10))
-XB_VAR_BV(:,:,:) = 0.
-ALLOCATE(XB_VAR_NOBV(NNB_TOPD_STEP,NNCAT,10))
-XB_VAR_NOBV(:,:,:) = 0.
-! 
-ALLOCATE(XB_VAR_TOT(NNB_TOPD_STEP,10))
-XB_VAR_TOT(:,:) = 0.
-!
-ALLOCATE(XB_RUNOFF_TOPD(KNI))
-ALLOCATE(XB_RUNOFF_TOPDM(KNI))
-ALLOCATE(XB_ATOP_BYMESH(KNI))
-!  
- CALL UNPACK_SAME_RANK(NR_NATURE,XATOP,XB_ATOP_BYMESH)
-IF(CRUNOFF=='TOPD')THEN
-  CALL UNPACK_SAME_RANK(NR_NATURE,XRUNOFF_TOP,XB_RUNOFF_TOPDM)
-ELSE
-  XB_RUNOFF_TOPDM = XB_RUNOFF_ISBAM
-ENDIF
-!
-ALLOCATE(YB_VARQ(5)) 
-YB_VARQ(1)='Q_TOT '
-YB_VARQ(2)='Q_RUN '
-YB_VARQ(3)='Q_DR  '
-YB_VARQ(4)='ST_RUN'
-YB_VARQ(5)='ST_DR '
-!
-!
-ALLOCATE(XB_QTOT(NNCAT))
-ALLOCATE(XB_QDR(NNCAT))
-ALLOCATE(XB_QRUN(NNCAT))
-ALLOCATE(XB_VAR_Q(NNB_TOPD_STEP,NNCAT,5))
-!
-ALLOCATE(XB_QTOTM(NNCAT))
-ALLOCATE(XB_QDRM(NNCAT))
-ALLOCATE(XB_QRUNM(NNCAT))
-! 
-!init var bilan q
-!
-XB_QTOT(:)=0.
-XB_QDR(:) =0.
-XB_QRUN(:)=0.
-!
-XB_VAR_Q(:,:,:)=0
-!
-DO JJ=1,NNCAT
-  XB_QTOTM(JJ) = SUM(XQTOT(JJ,:))
-  XB_QRUNM(JJ) = SUM(XQB_RUN(JJ,:))
-  XB_QDRM(JJ)  = SUM(XQB_DR(JJ,:))
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('INIT_BUDGET_COUPL_ROUT',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE INIT_BUDGET_COUPL_ROUT
diff --git a/src/SURFEX/init_coupl_topd.F90 b/src/SURFEX/init_coupl_topd.F90
deleted file mode 100644
index 2facfb7283e03a039902c7ef4de3e77918f7c172..0000000000000000000000000000000000000000
--- a/src/SURFEX/init_coupl_topd.F90
+++ /dev/null
@@ -1,471 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE INIT_COUPL_TOPD(HPROGRAM,KI)
-!     #######################
-!
-!!****  *INIT_COUPL_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!!     This routine aims at initialising the variables 
-!     needed for coupling with Topmodel.
-! 
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   16/10/2003
-!!      Modif BV : supression of variables specific to Topmodel
-!!      20/12/2007 - mll : Adaptation between a lonlat grid system for ISBA
-!!                         and lambert II projection for topmodel
-!!      11/2011: Modif BV : Creation of masks between ISBA and TOPODYN
-!                transfered in PGD step (routine init_pgd_topd)
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-! Modules
-USE MODD_COUPLING_TOPD, ONLY : XWSTOPI, XWFCTOPI, XDTOPI, XAS_NATURE, XATOP,&
-                                 XCSTOPI, XWTOPT, XAVG_RUNOFFCM, XAVG_DRAINCM,&
-                                 XDTOPT, XKA_PRE, XKAC_PRE, NMASKI, XDMAXFC, &
-                                 XWG_FULL, XWSTOPT, XWFCTOPT, NMASKT, & 
-                                 NNBV_IN_MESH, XBV_IN_MESH, XTOTBV_IN_MESH,&
-                                 XRUNOFF_TOP, NNPIX,&
-                                 XFRAC_D2, XFRAC_D3, XWGI_FULL,&
-                                 XRUN_TOROUT, XDR_TOROUT,&
-                                 LSTOCK_TOPD,NNB_STP_RESTART 
-USE MODD_DUMMY_EXP_PROFILE,ONLY :XF_PARAM, XC_DEPTH_RATIO
-USE MODD_TOPODYN,       ONLY : NNCAT, XMPARA, XCSTOPT, NMESHT, XDXT,&
-                                 NNMC, XRTOP_D2, NNB_TOPD_STEP
-!
-USE MODD_SURF_PAR,         ONLY : XUNDEF, NUNDEF
-USE MODD_ISBA_n,           ONLY : XSAND, XDG, XCLAY, XWG,&
-                                  CKSAT, XCONDSAT,XWGI, XF_PARAM_i=>XF_PARAM, &
-                                  XC_DEPTH_RATIO_i=>XC_DEPTH_RATIO
-USE MODD_DIAG_EVAP_ISBA_n, ONLY : XAVG_RUNOFFC, XAVG_DRAINC
-
-USE MODD_SURF_ATM_GRID_N,  ONLY : XMESH_SIZE
-USE MODD_SURF_ATM_n,       ONLY : NSIZE_NATURE, NR_NATURE
-!
-! Interfaces
-USE MODI_GET_LUOUT
-USE MODI_READ_FILE_MASKTOPD
-USE MODI_PACK_SAME_RANK
-USE MODI_UNPACK_SAME_RANK
-USE MODI_ISBA_TO_TOPD
-USE MODI_RESTART_COUPL_TOPD
-!
-USE MODE_SOIL
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*), INTENT(IN) :: HPROGRAM    ! 
-INTEGER, INTENT(IN)          :: KI          ! Grid dimensions
-!
-!*      0.2    declarations of local variables
-!
-!REAL, DIMENSION(:), ALLOCATABLE   :: ZDTAV                            ! Averaged depth soil on TOP-LAT grid
-REAL, DIMENSION(:), ALLOCATABLE   :: ZSAND_FULL, ZCLAY_FULL, ZDG_FULL ! Isba variables on the full domain
-REAL, DIMENSION(:), ALLOCATABLE   :: ZFRAC    ! fraction of SurfEx mesh that covers one or several catchments
-REAL, DIMENSION(:), ALLOCATABLE   :: ZDMAXAV  ! dificit maximal moyen par bassin
-REAL, DIMENSION(:),ALLOCATABLE    :: ZSANDTOPI, ZCLAYTOPI!, ZWWILTTOPI !sand and clay fractions on TOPMODEL layers
-!
-!ludo
-REAL, DIMENSION(:), ALLOCATABLE   :: ZKSAT       !ksat surf 
-REAL, DIMENSION(:), ALLOCATABLE   :: ZF_PARAM_FULL
-REAL, DIMENSION(:,:), ALLOCATABLE :: ZF_PARAMT!, ZWWILTTOPT
-REAL, DIMENSION(:), ALLOCATABLE   :: ZDG2_FULL, ZDG3_FULL, ZWG2_FULL, ZWG3_FULL, ZRTOP_D2
-REAL,DIMENSION(:), ALLOCATABLE    :: ZWGI_FULL, Z_WFCTOPI, Z_WSTOPI
-!
-REAL                              :: ZCOEF_ANIS  !coefficient anisotropie Ksat:
-                                                 ! Ksat horiz=ZCOEF*Ksat vert
-INTEGER                   :: JJ,JI           ! loop control 
-INTEGER                   :: JCAT,JMESH      ! loop control 
-INTEGER                   :: ILUOUT          ! Logical unit for output filr
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('INIT_COUPL_TOPD',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-WRITE(ILUOUT,*) 'INITIALISATION INIT_COUPL_TOPD'
-!
-ALLOCATE(NMASKT(NNCAT,NMESHT))
-NMASKT(:,:) = NUNDEF
-!
-!*       1    Initialization:
-!               ---------------
-!
-! la surface saturee, ā l'initialisation est nulle, donc on initialise les lambdas de telle sorte qu'aucun pixel ne soit sature
-ALLOCATE(XKA_PRE (NNCAT,NMESHT))
-ALLOCATE(XKAC_PRE(NNCAT))
-XKA_PRE(:,:) = 0.0
-XKAC_PRE(:)  = MAXVAL(XKA_PRE) + 1.
-!
-!Cumulated runoff initialisation
-ALLOCATE(XRUNOFF_TOP(NSIZE_NATURE))
-XRUNOFF_TOP  (:) = XAVG_RUNOFFC(:)
-!
-IF(.NOT.ALLOCATED(XAVG_RUNOFFCM)) ALLOCATE(XAVG_RUNOFFCM(NSIZE_NATURE))
-XAVG_RUNOFFCM(:) = XAVG_RUNOFFC(:)
-!
-IF(.NOT.ALLOCATED(XAVG_DRAINCM )) ALLOCATE(XAVG_DRAINCM (NSIZE_NATURE))
-XAVG_DRAINCM (:) = XAVG_DRAINC(:)
-!
-!
-! Reading masks
- CALL READ_FILE_MASKTOPD(KI)
-!
-!*      2.1     Fraction of SurfEx mesh with TOPMODEL
-!               -------------------------------------
-!
-ALLOCATE(NNBV_IN_MESH  (KI,NNCAT))
-ALLOCATE(XBV_IN_MESH   (KI,NNCAT))
-ALLOCATE(XTOTBV_IN_MESH(KI))
-!
-XTOTBV_IN_MESH(:) = 0.0
-!
-DO JJ=1,KI
-  !
-  XBV_IN_MESH(JJ,:) = 0.0
-  !
-  DO JI=1,NNCAT
-    NNBV_IN_MESH(JJ,JI) = COUNT( NMASKI(JJ,JI,:)/=NUNDEF )
-    XBV_IN_MESH (JJ,JI) = REAL(NNBV_IN_MESH(JJ,JI)) * XDXT(JI)**2
-    XTOTBV_IN_MESH (JJ) = XTOTBV_IN_MESH(JJ) + XBV_IN_MESH(JJ,JI)
-  ENDDO
-  !
-ENDDO
-!
-!*      2.2     Fraction of SurfEx mesh with each catchment
-!               -------------------------------------------
-!
-ALLOCATE(ZFRAC(KI))  ! fraction not covered by catchments
-ZFRAC(:) = ( XMESH_SIZE(:)-XTOTBV_IN_MESH(:) ) / XMESH_SIZE(:)
-ZFRAC(:) = MIN(MAX(ZFRAC(:),0.),1.)
-!
-ALLOCATE(XATOP(NSIZE_NATURE)) ! fraction covered by catchments part nature
- CALL PACK_SAME_RANK(NR_NATURE,(1.-ZFRAC),XATOP)
-!
-!
-IF (HPROGRAM=='POST  ') GOTO 10
-!
-!*      3.0     Wsat, Wfc and depth for TOPODYN on ISBA grid
-!               --------------------------------------------
-!*      3.1     clay, sand fraction, depth hydraulic conductivity at saturation of the layer for TOPODYN
-!               ---------------------------------------------------------
-!
-ALLOCATE(ZSAND_FULL(KI))
-ALLOCATE(ZCLAY_FULL(KI))
- CALL UNPACK_SAME_RANK(NR_NATURE,XSAND(:,2),ZSAND_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XCLAY(:,2),ZCLAY_FULL)
-!
-!ludo prof variable pour tr lat (OK car sol homogene verticalement, faux sinon)
-ALLOCATE(ZDG2_FULL(KI))
-ALLOCATE(ZDG3_FULL(KI))
- CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,2,1),ZDG2_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XDG(:,3,1),ZDG3_FULL)
-!
-ALLOCATE(ZRTOP_D2(KI))
-ZRTOP_D2(:) = 0.
-!
-DO JMESH=1,KI
-  IF ( ZDG2_FULL(JMESH)/=XUNDEF .AND. ZFRAC(JMESH)<1. ) THEN
-    ZRTOP_D2(JMESH) = 0.
-    DO JCAT=1,NNCAT
-     !moyenne ponderee pour cas ou plusieurs BV sur maille
-      ZRTOP_D2(JMESH) = ZRTOP_D2(JMESH) + XRTOP_D2(JCAT)*XBV_IN_MESH(JMESH,JCAT)/XMESH_SIZE(JMESH)    
-    END DO
-  ENDIF   
-ENDDO
-!ZTOP_D2 * D2 < D3 : the depth concernet by lateral transfers is lower than D2
-WHERE( ZDG2_FULL/=XUNDEF .AND. ZRTOP_D2*ZDG2_FULL>ZDG3_FULL ) ZRTOP_D2(:) = ZDG3_FULL(:)/ZDG2_FULL(:)
-!
-DEALLOCATE(ZFRAC)
-!
-ALLOCATE(XFRAC_D2 (KI))
-ALLOCATE(XFRAC_D3 (KI))
-XFRAC_D2(:)=1.
-XFRAC_D3(:)=0.
-WHERE( ZDG2_FULL/=XUNDEF .AND. ZRTOP_D2*ZDG2_FULL>ZDG2_FULL  ) ! if the depth is > D2
-  XFRAC_D2(:) = MIN(1.,ZRTOP_D2(:))
-  XFRAC_D3(:) = (ZRTOP_D2(:)*ZDG2_FULL(:)-ZDG2_FULL(:)) / (ZDG3_FULL(:)-ZDG2_FULL(:))
-  XFRAC_D3(:) = MAX(0.,XFRAC_D3(:))
-END WHERE
-!
-ALLOCATE(ZDG_FULL(KI))
-WHERE (ZDG2_FULL/=XUNDEF)
-  ZDG_FULL = XFRAC_D2*ZDG2_FULL + XFRAC_D3*(ZDG3_FULL-ZDG2_FULL)
-ELSEWHERE
-  ZDG_FULL = XUNDEF
-END WHERE
-!
-ALLOCATE(ZSANDTOPI(KI))
-ALLOCATE(ZCLAYTOPI(KI))
-ZSANDTOPI(:)=0.0
-ZCLAYTOPI(:)=0.0
-ALLOCATE(XDTOPI(KI))
-XDTOPI(:)=0.0
-WHERE ( ZDG_FULL/=XUNDEF .AND. ZDG_FULL/=0. )
-  XDTOPI = ZDG_FULL
-  ZSANDTOPI = ZSANDTOPI + ZSAND_FULL * ZDG_FULL
-  ZCLAYTOPI = ZCLAYTOPI + ZCLAY_FULL * ZDG_FULL
-  ZSANDTOPI = ZSANDTOPI / XDTOPI
-  ZCLAYTOPI = ZCLAYTOPI / XDTOPI
-ELSEWHERE
-  ZSANDTOPI = XUNDEF
-  ZCLAYTOPI = XUNDEF
-  XDTOPI = XUNDEF
-END WHERE
-DEALLOCATE(ZSAND_FULL)
-DEALLOCATE(ZCLAY_FULL)
-!
-!*      4.1     depth of the Isba layer on TOP-LAT grid
-!               ---------------------------------------
-!
-ALLOCATE(XDTOPT(NNCAT,NMESHT))
-XDTOPT(:,:) = 0.0
- CALL ISBA_TO_TOPD(XDTOPI,XDTOPT)
-!
-!*      3.2     Wsat and Wfc on TOPODYN layer
-!               -----------------------------
-!
-ALLOCATE(XWSTOPI   (KI))
-ALLOCATE(XWFCTOPI  (KI))
-XWSTOPI (:) = 0.0
-XWFCTOPI(:) = 0.0
-!ALLOCATE(ZWWILTTOPI(KI))
-XWSTOPI    = WSAT_FUNC_1D (ZCLAYTOPI,ZSANDTOPI,'CH78')
-XWFCTOPI   = WFC_FUNC_1D  (ZCLAYTOPI,ZSANDTOPI,'CH78')
-!ZWWILTTOPI = WWILT_FUNC_1D(ZCLAYTOPI,ZSANDTOPI,'CH78')
-!
-!modif ludo test ksat exp
-WRITE(ILUOUT,*) 'CKSAT==',CKSAT
-
-ALLOCATE(ZKSAT(KI))
-ZKSAT   (:) = 0.0
-ALLOCATE(XCSTOPI(KI))
-XCSTOPI(:) = 0.0
-IF( CKSAT=='SGH' .OR. CKSAT=='EXP' ) THEN
-  !
-  !ludo calcul des profondeur efficaces
-  !ZRTOP_D2(:) = 1.
-  !ALLOCATE(XC_DEPTH_RATIO(SIZE(XC_DEPTH_RATIO_i)))
-  !XC_DEPTH_RATIO(:) = XC_DEPTH_RATIO_i(:)
-  !ZRTOP_D2(:) = XC_DEPTH_RATIO(:)
-  !valeur patch 1 (idem wsat wfc) a voir cas ou il existe plusieurs patchs 
-  CALL UNPACK_SAME_RANK(NR_NATURE,XCONDSAT(:,1,1),ZKSAT)
-  !passage de definition Ksat(profondeur) en Ksat(deficit)
-  WHERE ( ZDG_FULL/=XUNDEF .AND. (XWSTOPI-XWFCTOPI/=0.) )
-    XCSTOPI(:) = ZKSAT(:) / (XWSTOPI(:)-XWFCTOPI(:))
-  END WHERE
-  !
-ELSE
-  !
-  XCSTOPI(:) = HYDCONDSAT_FUNC(ZCLAYTOPI,ZSANDTOPI,'CH78')
-  !passage de definition Ksat(profondeur) en Ksat(deficit)
-  WHERE ( ZDG_FULL/=XUNDEF .AND. (XWSTOPI-XWFCTOPI/=0.) )
-    XCSTOPI(:) = XCSTOPI(:) / (XWSTOPI(:)-XWFCTOPI(:))
-  END WHERE
-  !
-ENDIF
-!
-DEALLOCATE(ZSANDTOPI)
-DEALLOCATE(ZCLAYTOPI)
-DEALLOCATE(ZRTOP_D2)
-!
-!*      4.3     Ko on TOP-LAT grid
-!               ------------------
-!
-ALLOCATE(XCSTOPT(NNCAT,NMESHT))
- CALL ISBA_TO_TOPD(XCSTOPI,XCSTOPT)
-WHERE (XCSTOPT == XUNDEF) XCSTOPT = 0.0
-!
-!*      3.3     Initialization of the previous time step water storage on ISBA grid to calculate the refill on Isba grid
-!               -------------------------------------------------------------------------
-!
-ALLOCATE(ZWG2_FULL(KI))
-ALLOCATE(ZWG3_FULL(KI))
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,2,1),ZWG2_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XWG(:,3,1),ZWG3_FULL)
-!
-ALLOCATE(XWG_FULL(KI))
-WHERE ( ZWG2_FULL/=XUNDEF .AND. ZDG_FULL/=0. )
-  XWG_FULL = ( XFRAC_D2*ZDG2_FULL*ZWG2_FULL + XFRAC_D3*(ZDG3_FULL-ZDG2_FULL)*ZWG3_FULL ) / ZDG_FULL
-ELSEWHERE
-  XWG_FULL = XUNDEF
-END WHERE
-!
-DEALLOCATE(ZWG2_FULL)
-DEALLOCATE(ZWG3_FULL)
-DEALLOCATE(ZDG3_FULL)
-!
-!
-ALLOCATE(XWTOPT(NNCAT,NMESHT))
-XWTOPT(:,:) = 0.0
- CALL ISBA_TO_TOPD(XWG_FULL,XWTOPT)
-WHERE (XWTOPT == XUNDEF) XWTOPT = 0.0
-!
-!ludo prise en compte glace (pas de glace dans 3e couche)
-ALLOCATE(ZWGI_FULL(KI))
-ALLOCATE(XWGI_FULL(KI))
- CALL UNPACK_SAME_RANK(NR_NATURE,XWGI(:,2,1),ZWGI_FULL)
-!
-WHERE ( ZWGI_FULL/=XUNDEF .AND. XFRAC_D2>0. .AND. ZDG_FULL/=0. )
-  XWGI_FULL = (XFRAC_D2*ZDG2_FULL*ZWGI_FULL) / ZDG_FULL
-ELSEWHERE
-  XWGI_FULL = XUNDEF
-END WHERE
-!
-DEALLOCATE(ZWGI_FULL)
-DEALLOCATE(ZDG2_FULL)
-DEALLOCATE(ZDG_FULL)
-!
-ALLOCATE(Z_WFCTOPI(KI))
-ALLOCATE(Z_WSTOPI (KI))
-!test reservoir top=eau+glace -> pas de modif Wsat et Wfc
-WHERE ( XWGI_FULL/=XUNDEF .AND. XWSTOPI/=0. )
-  Z_WSTOPI  = XWSTOPI - XWGI_FULL
-  Z_WFCTOPI = XWFCTOPI * Z_WSTOPI / XWSTOPI
-END WHERE
-!
-!ludo calcul en fct teneur glace
-!
-ALLOCATE(XWSTOPT (NNCAT,NMESHT))
-ALLOCATE(XWFCTOPT(NNCAT,NMESHT))
- CALL ISBA_TO_TOPD(Z_WSTOPI,XWSTOPT)
- CALL ISBA_TO_TOPD(Z_WFCTOPI,XWFCTOPT)
-DEALLOCATE(Z_WSTOPI)
-DEALLOCATE(Z_WFCTOPI)
-!
-!*      4.0     calcul of time constant variables on TOPODYN grid 
-!               -------------------------------------------------
-!
-!*      4.2     Wsat and Wfc on TOP-LAT grid
-!               ----------------------------
-!
-ALLOCATE(XDMAXFC(NNCAT,NMESHT))
-XDMAXFC(:,:) = XUNDEF
-WHERE (XWFCTOPT /= XUNDEF) XDMAXFC = (XWSTOPT - XWFCTOPT) * XDTOPT ! (m)
-!
-!
-!*      4.4     Initialisation of the previous time step water storage on topodyn-lat grid
-!               --------------------------------------------------------------------------
-!*      4.5     M parameter on TOPODYN grid
-!               ------------------------
-!*      4.5.1   Mean depth soil on catchment
-!
-ALLOCATE(XMPARA (NNCAT))
-XMPARA  (:) = 0.0
-!
-IF( CKSAT=='EXP' .OR. CKSAT=='SGH' ) THEN
-  !ludo test
-  ALLOCATE(ZF_PARAM_FULL(KI))
-  ALLOCATE(ZF_PARAMT(NNCAT,NMESHT))
-  ALLOCATE(XF_PARAM(SIZE(XF_PARAM_i)))
-  XF_PARAM(:) = XF_PARAM_i(:)
-  CALL UNPACK_SAME_RANK(NR_NATURE,XF_PARAM(:),ZF_PARAM_FULL)
-  CALL ISBA_TO_TOPD(ZF_PARAM_FULL,ZF_PARAMT)
-  DEALLOCATE(ZF_PARAM_FULL)
-  !
-  !passage de f a M (M=Wsat-Wfc/f)
-  !ludo test ksat exp
-  !ALLOCATE(ZWWILTTOPT(NNCAT,NMESHT))
-  !CALL ISBA_TO_TOPD(ZWWILTTOPI,ZWWILTTOPT)
-  WHERE( ZF_PARAMT/=XUNDEF .AND. ZF_PARAMT/=0. ) ZF_PARAMT = (XWSTOPT-XWFCTOPT)/ZF_PARAMT
-  !DEALLOCATE(ZWWILTTOPT)
-  !
-  DO JJ=1,NNCAT
-    XMPARA(JJ) = SUM(ZF_PARAMT(JJ,:),MASK=ZF_PARAMT(JJ,:)/=XUNDEF) / NNMC(JJ)
-  ENDDO
-  !
-  ZCOEF_ANIS = 1.
-  XCSTOPT = XCSTOPT*ZCOEF_ANIS
-  !
-  DEALLOCATE(ZF_PARAMT)
-  !
-ELSE
-  !
-  ALLOCATE(ZDMAXAV(NNCAT))
-  ZDMAXAV(:) = 0.0
-  DO JJ=1,NNCAT
-    ZDMAXAV(JJ) = SUM( XDMAXFC(JJ,:),MASK=XDMAXFC(JJ,:)/=XUNDEF ) / NNMC(JJ)
-  ENDDO
-  !
-  !ALLOCATE(ZDTAV  (NNCAT))
-  !ZDTAV   (:) = 0.0
-  DO JJ=1,NNCAT 
-    !ZDTAV(JJ) = SUM(XDTOPT(JJ,:),MASK=XDTOPT(JJ,:)/=XUNDEF) / NNMC(JJ)
-    XMPARA(JJ) = ZDMAXAV(JJ) / 4.
-  ENDDO
-  !DEALLOCATE(ZDTAV)
-  DEALLOCATE(ZDMAXAV)
-  !
-ENDIF
-!
-!DEALLOCATE(ZWWILTTOPI)
-! 
-!*      5.0      Initial saturated area computation
-!               -----------------------------------------------------------
-!
-ALLOCATE(XAS_NATURE(NSIZE_NATURE))
-XAS_NATURE(:) = 0.0
-!
-!*      6.0     Stock management in case of restart
-!               -----------------------------------------------------------
-!
-10 CONTINUE
-!
-!stock
-ALLOCATE(XRUN_TOROUT(NNCAT,NNB_TOPD_STEP+NNB_STP_RESTART))
-ALLOCATE(XDR_TOROUT (NNCAT,NNB_TOPD_STEP+NNB_STP_RESTART))
-XRUN_TOROUT(:,:) = 0.
-XDR_TOROUT (:,:) = 0.
-!
-IF (HPROGRAM=='POST  ') GOTO 20
-!
-IF (LSTOCK_TOPD) CALL RESTART_COUPL_TOPD(HPROGRAM,KI)
-!
-!*      7.0     deallocate
-!               ----------
-!
-20 CONTINUE
-!
-IF (LHOOK) CALL DR_HOOK('INIT_COUPL_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE INIT_COUPL_TOPD
-
-
-
-
-
-
-
diff --git a/src/SURFEX/init_surf_topd.F90 b/src/SURFEX/init_surf_topd.F90
deleted file mode 100644
index 3d44a35953a0759bcc12eb4ba7c8754280e74ed4..0000000000000000000000000000000000000000
--- a/src/SURFEX/init_surf_topd.F90
+++ /dev/null
@@ -1,113 +0,0 @@
-!-------------------------------------------------------------------------------
-!     #############################################################
-      SUBROUTINE INIT_SURF_TOPD(HPROGRAM,KI)
-!     #############################################################
-!
-!!****  *INIT_SURF_TOPD* - routine to initialize variables needed for coupling with Topmodel
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!**  METHOD
-!!    ------
-!!    The routine open and read the namelists NAM_COUPL_TOPD and NAM_TOPD,
-!! calculates the number of catchments concerned, the different time step 
-!! variables and all the variables nedded for coupling with Topmodel. 
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!	B. Vincendon   *Meteo France*	
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    11/2006
-!!      Modif 04/2007: Arguments PTOPD_STEP,KNB_TOPD_STEP become module
-!!                     variables from MODD_TOPDDYN
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-USE MODD_SURFEX_MPI, ONLY : NPROC
-USE MODD_SURFEX_OMP, ONLY : NBLOCKTOT
-!
-USE MODD_TOPODYN, ONLY :CCAT, XSPEEDR, XSPEEDH, NNCAT, &
-                        XRTOP_D2, XSPEEDG
-USE MODD_COUPLING_TOPD, ONLY :  LCOUPL_TOPD, NNB_TOPD, LBUDGET_TOPD
-!
-USE MODD_ISBA_n, ONLY : TSNOW
-!
-USE MODI_GET_LUOUT
-USE MODI_ABOR1_SFX
-!
-USE MODI_INIT_TOPD
-USE MODI_INIT_COUPL_TOPD
-USE MODI_INIT_BUDGET_COUPL_ROUT
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of arguments
-!              -------------------------
-!
- CHARACTER(LEN=*),  INTENT(IN)     :: HPROGRAM      !
-INTEGER,           INTENT(IN)     :: KI            ! grid dimension
-!
-INTEGER :: ILUOUT
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('INIT_SURF_TOPD',0,ZHOOK_HANDLE)
-!
-IF (LCOUPL_TOPD) THEN
-  IF (NPROC>1) CALL ABOR1_SFX('INIT_SURF_TOPD: TOPD CANNOT RUN WITH MORE THAN 1 MPI TASK') 
-  IF (NBLOCKTOT>1) CALL ABOR1_SFX("INIT_SURF_TOPD: TOPD CANNOT RUN WITH NUMEROUS OPENMP BLOCKS")
-  IF (TSNOW%SCHEME/='3-L') &
-        CALL ABOR1_SFX("INIT_SURF_TOPD: coupling with topmodel only runs with TSNOW%SCHEME=3-L")
-ENDIF
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!            
-!         1.   Reads the namelists
-!              --------------------
-!
-WRITE(ILUOUT,*) 'Debut init_surf_topo_n'
-!
-IF (LCOUPL_TOPD) THEN
-  !
-  !         3.   Initialises variables specific to Topmodel
-  !              -------------------------------------------
-  WRITE(ILUOUT,*) 'NNCAT',NNCAT
-  !
-  CALL INIT_TOPD(HPROGRAM)
-  !
-  !         4.   Initialises variables nedded for coupling with Topmodel
-  !              -------------------------------------------------------
-  !
-  CALL INIT_COUPL_TOPD(HPROGRAM,KI)
-  !
-  WRITE(ILUOUT,*) 'Couplage avec TOPMODEL active'
-  !
-  IF (LBUDGET_TOPD) CALL INIT_BUDGET_COUPL_ROUT(KI)
-  !
-ELSE
-  !
-  WRITE(ILUOUT,*) 'Pas de couplage avec TOPMODEL'
-  !
-ENDIF
-!
-IF (LHOOK) CALL DR_HOOK('INIT_SURF_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE INIT_SURF_TOPD
diff --git a/src/SURFEX/init_topd.F90 b/src/SURFEX/init_topd.F90
deleted file mode 100644
index 51b3e23a11dabf49029594e016ddd6735ef35532..0000000000000000000000000000000000000000
--- a/src/SURFEX/init_topd.F90
+++ /dev/null
@@ -1,276 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE INIT_TOPD(HPROGRAM)
-!     #######################
-!
-!!****  *INIT_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!     This routine aims at initialising the variables 
-!     needed of running Topmodel.
-!              
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR 
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!      B. Vincendon    * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2006
-!!      Modification 04/2007 : Supression of 2 arguments KSTEP,PSTEP 
-!!                             that are now module arguments from MODD_TOPDDYN_n
-!!                             NNB_TOPD_STEP,XTOPD_STEP
-!!      Modification 11/2011 : Exfiltration option removed (B. Vincendon)
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_COUPLING_TOPD, ONLY : NNB_STP_RESTART
-USE MODD_TOPODYN,       ONLY : CCAT, NNCAT, NNB_TOPD_STEP, XTOPD_STEP,&
-                               XDXT, NNXC, NNYC,&
-                               XNUL, XX0, XY0, NNPT,&
-                               NX_STEP_ROUT, XSPEEDR,&
-                               XSPEEDH, NNMC, NMESHT, NPMAX,&
-                               NLINE,  XDMAXT,&
-                               XTOPD, XDRIV, XDHIL, XTIME_TOPD,&
-                               XDGRD, XSPEEDG, XTIME_TOPD_DRAIN,&
-                               XQTOT, XTANB, XSLOP, XDAREA,&
-                               XLAMBDA, XCONN, XQB_DR, XQB_RUN
-!
-USE MODD_TOPD_PAR, ONLY : NDIM
-USE MODD_SURF_PAR,       ONLY : XUNDEF, NUNDEF
-!
-USE MODI_GET_LUOUT
-USE MODI_READ_TOPD_HEADER_DTM
-USE MODI_READ_TOPD_FILE
-USE MODI_READ_TOPD_HEADER_CONNEX
-USE MODI_READ_CONNEX_FILE
-USE MODI_READ_SLOPE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*), INTENT(IN) :: HPROGRAM    !
-!
-!*      0.2    declarations of local variables
-!
-!
- CHARACTER(LEN=50), DIMENSION(:),ALLOCATABLE :: YFILETOP ! topographic file names
- CHARACTER(LEN=50), DIMENSION(:),ALLOCATABLE :: YFILECON ! topographic file names
- CHARACTER(LEN=50), DIMENSION(:),ALLOCATABLE :: YFILESLO ! topographic file names
- CHARACTER(LEN=50), DIMENSION(:),ALLOCATABLE :: YFILEDH  ! topographic file names
- CHARACTER(LEN=50), DIMENSION(:),ALLOCATABLE :: YFILEDR  ! topographic file names
-INTEGER                   :: JJ,JCAT ! loop control 
-INTEGER                   :: IUNIT                  ! Unit of the files
-INTEGER                   :: ILUOUT                 ! Unit of the files
-!
-REAL, DIMENSION(:),ALLOCATABLE    :: ZTOPD_READ !Topgraphic variable read
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('INIT_TOPD',0,ZHOOK_HANDLE)
-!
-!*       1    Initialization:
-!               ---------------
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-WRITE(ILUOUT,*) 'INITIALISATION INIT_TOPD'
-!
-ALLOCATE(YFILETOP(NNCAT))
-ALLOCATE(YFILECON(NNCAT))
-ALLOCATE(YFILESLO(NNCAT))
-ALLOCATE(YFILEDH (NNCAT))
-ALLOCATE(YFILEDR (NNCAT))
-!
-DO JCAT=1,NNCAT
-  YFILETOP(JCAT)=TRIM(CCAT(JCAT))//'_FilledDTM.map'
-  YFILECON(JCAT)=TRIM(CCAT(JCAT))//'_connections.vec'
-  YFILESLO(JCAT)=TRIM(CCAT(JCAT))//'_slope.vec'
-  YFILEDR(JCAT)=TRIM(CCAT(JCAT))//'_RiverDist.map'
-  YFILEDH(JCAT)=TRIM(CCAT(JCAT))//'_HillDist.map'
-ENDDO
-
-ALLOCATE(NNMC(NNCAT))
-NNMC(:)=0
-ALLOCATE(NNXC(NNCAT))
-NNXC(:)=0
-ALLOCATE(NNYC(NNCAT))
-NNYC(:)=0
-ALLOCATE(XX0(NNCAT))
-XX0(:)=0.0
-ALLOCATE(XY0(NNCAT))
-XY0(:)=0.0
-ALLOCATE(NNPT(NNCAT))
-NNPT(:)=0
-ALLOCATE(XNUL(NNCAT))
-ALLOCATE(XDXT(NNCAT))
-!
-!*       2      Topographic files
-!               -----------------------------
-DO JCAT=1,NNCAT
-  CALL READ_TOPD_HEADER_DTM(HPROGRAM,YFILETOP(JCAT),'FORMATTED',&
-                            XX0(JCAT),XY0(JCAT),NNXC(JCAT),NNYC(JCAT),&
-                            XNUL(JCAT),XDXT(JCAT))
-ENDDO
-!
-NNPT(:) = NNXC(:) * NNYC(:)
-!
-NPMAX = MAXVAL(NNPT(:))
-!
-!
-ALLOCATE(NLINE(NNCAT,NPMAX))
-NLINE(:,:)=0
-ALLOCATE(XDRIV(NNCAT,NPMAX))
-XDRIV(:,:)=0.0
-ALLOCATE(XDHIL(NNCAT,NPMAX))
-XDHIL(:,:)=0.0
-ALLOCATE(XDGRD(NNCAT,NPMAX))
-XDGRD(:,:)=0.0
-ALLOCATE(XQTOT(NNCAT,NNB_TOPD_STEP))
-XQTOT(:,:)=0.0
-ALLOCATE(XTOPD(NNCAT,NPMAX))
-XTOPD(:,:)=0.0
-
-ALLOCATE(ZTOPD_READ(NPMAX))
-DO JCAT = 1,NNCAT
-  !
-  CALL READ_TOPD_FILE(HPROGRAM,YFILETOP(JCAT),'FORMATTED',NNPT(JCAT),ZTOPD_READ)
-  DO JJ = 1,NNPT(JCAT)
-    XTOPD(JCAT,JJ) = ZTOPD_READ(JJ) ! XTOPD can only be >=0
-  ENDDO
-  !
-ENDDO
-! for XTOPD, we do not have to use NLINE, all the pixels of the rectancle around
-! the catchment are read.
-!
-!*      4       Connection files
-!               ----------------
-!
-DO JCAT=1,NNCAT
-  CALL READ_TOPD_HEADER_CONNEX(HPROGRAM,YFILECON(JCAT),'FORMATTED',NNMC(JCAT))
-ENDDO
-!
-NMESHT=MAXVAL(NNMC(:))
-!
-ALLOCATE(XDMAXT(NNCAT,NMESHT))
-ALLOCATE(XCONN(NNCAT,NMESHT,NDIM))
-XCONN(:,:,:)=0.0
-XDMAXT(:,:)=0.0
-!
-DO JCAT=1,NNCAT
-  CALL READ_CONNEX_FILE(HPROGRAM,YFILECON(JCAT),'FORMATTED',NNMC(JCAT),XCONN(JCAT,:,:),NLINE(JCAT,:))
-ENDDO
-!
-!*     5        Slope files
-!               -----------
-IF (HPROGRAM/='PGD') THEN
-  !
-  ALLOCATE(XTANB(NNCAT,NMESHT))
-  ALLOCATE(XSLOP(NNCAT,NMESHT))
-  ALLOCATE(XDAREA(NNCAT,NMESHT))
-  ALLOCATE(XLAMBDA(NNCAT,NMESHT))
-  !
-  DO JCAT=1,NNCAT
-    CALL READ_SLOPE_FILE(HPROGRAM,YFILESLO(JCAT),'FORMATTED',NNMC(JCAT),&
-                       XTANB(JCAT,:),XSLOP(JCAT,:),XDAREA(JCAT,:),XLAMBDA(JCAT,:))
-  ENDDO
-  !
-  !*      6       River Distance file
-  !               -------------------
-  !
-  DO JCAT=1,NNCAT
-    CALL  READ_TOPD_FILE(HPROGRAM,YFILEDR(JCAT),'FORMATTED',NNPT(JCAT),ZTOPD_READ)
-    DO JJ=1,NPMAX
-      IF ( NLINE(JCAT,JJ)/=0. .AND. NLINE(JCAT,JJ)/=XUNDEF ) &
-       XDRIV(JCAT,NLINE(JCAT,JJ)) = ZTOPD_READ(JJ)
-    ENDDO
-  ENDDO
-  ! for XDRIV, we must use NLINE, online pixels inside de the catchment are read.
-  !
-  !*      7       Hillslope Distance file
-  !  -----------------------
-  !
-  DO JCAT=1,NNCAT
-    CALL  READ_TOPD_FILE(HPROGRAM,YFILEDH(JCAT),'FORMATTED',NNPT(JCAT),ZTOPD_READ)
-    DO JJ=1,NPMAX
-      IF ( NLINE(JCAT,JJ)/=0. .AND. NLINE(JCAT,JJ)/=XUNDEF ) &
-        XDHIL(JCAT,NLINE(JCAT,JJ)) = ZTOPD_READ(JJ)
-    ENDDO
-  ENDDO
-  !
-  ! for XDHIL, we must use NLINE, online pixels inside de the catchment are read.
-  XDGRD(:,:) = XDHIL(:,:)
-  !
-  !*      8       Calculations for routing by geomorpho
-  !               -------------------------------------
-  !
-  ALLOCATE(NX_STEP_ROUT(NNCAT))
-  ALLOCATE(XTIME_TOPD(NNCAT,NMESHT))
-  ALLOCATE(XTIME_TOPD_DRAIN(NNCAT,NMESHT))
-  !
-  XTIME_TOPD(:,:) = 0.0
-  XTIME_TOPD_DRAIN(:,:) = 0.0
-  XSPEEDH(:) = XSPEEDR(:)/10.
-  !
-  !
-  DO JCAT=1,NNCAT
-    !
-    IF ( XSPEEDR(JCAT)/=0. .AND. XSPEEDG(JCAT)/=0. ) THEN
-      WHERE ( XDHIL(JCAT,1:NNMC(JCAT))/=XUNDEF .AND. XDRIV(JCAT,1:NNMC(JCAT))/=XUNDEF ) &
-        XTIME_TOPD(JCAT,1:NNMC(JCAT)) = XDHIL(JCAT,1:NNMC(JCAT)) / XSPEEDH(JCAT) + &
-                                        XDRIV(JCAT,1:NNMC(JCAT)) / XSPEEDR(JCAT)
-      WHERE ( XDGRD(JCAT,1:NNMC(JCAT))/=XUNDEF .AND. XDRIV(JCAT,1:NNMC(JCAT))/=XUNDEF ) &
-        XTIME_TOPD_DRAIN(JCAT,1:NNMC(JCAT)) = XDGRD(JCAT,1:NNMC(JCAT)) / XSPEEDG(JCAT) + &
-                                              XDRIV(JCAT,1:NNMC(JCAT)) / XSPEEDR(JCAT)
-    ELSE 
-      WRITE(ILUOUT,*) 'You have to choose some values for routing velocities'
-    ENDIF
-    !
-    IF (XTOPD_STEP/=0.) &
-      NX_STEP_ROUT(JCAT) = INT(MAXVAL(XTIME_TOPD(JCAT,1:NNMC(JCAT))) / XTOPD_STEP) + 1
-    !
-  ENDDO
-  !
-  IF ( NNB_STP_RESTART==0 ) NNB_STP_RESTART = MAX(NNB_TOPD_STEP,MAXVAL(NX_STEP_ROUT(:)))
-  !
-  !
-  ALLOCATE(XQB_DR(NNCAT,NNB_TOPD_STEP))
-  XQB_DR(:,:)=0.0
-  ALLOCATE(XQB_RUN(NNCAT,NNB_TOPD_STEP))
-  XQB_RUN(:,:)=0.0
-  !
-ENDIF
-!
-IF (LHOOK) CALL DR_HOOK('INIT_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE INIT_TOPD
diff --git a/src/SURFEX/isba_to_topd.F90 b/src/SURFEX/isba_to_topd.F90
deleted file mode 100644
index a28e9ee9b238356d75f17fd44758f6fa411888cf..0000000000000000000000000000000000000000
--- a/src/SURFEX/isba_to_topd.F90
+++ /dev/null
@@ -1,76 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE ISBA_TO_TOPD(PVARI,PVART)
-!     #######################
-!
-!!****  *ISBA_TO_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!
-!     
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   12/2003
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, NNMC
-USE MODD_COUPLING_TOPD, ONLY : NMASKT
-USE MODD_SURF_PAR,        ONLY : XUNDEF, NUNDEF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-REAL, DIMENSION(:), INTENT(IN)      :: PVARI   ! variable from ISBA grid
-REAL, DIMENSION(:,:), INTENT(OUT)   :: PVART   ! variable for TOPODYN grid
-!
-!*      0.2    declarations of local variables
-!
-INTEGER            :: JJ, JI             ! loop control 
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('ISBA_TO_TOPD',0,ZHOOK_HANDLE)
-!
-!*       1.     ISBA => TOPODYN-LAT
-!               -------------------
-!
-PVART(:,:)=XUNDEF
-DO JJ=1,NNCAT
-  DO JI=1,NNMC(JJ)
-    IF (NMASKT(JJ,JI)/=NUNDEF) PVART(JJ,JI) = PVARI(NMASKT(JJ,JI))
-  ENDDO
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('ISBA_TO_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE ISBA_TO_TOPD
diff --git a/src/SURFEX/isba_to_topdsat.F90 b/src/SURFEX/isba_to_topdsat.F90
deleted file mode 100644
index b1f2d90667adfa19bf43a3cb9cecc63646f5b365..0000000000000000000000000000000000000000
--- a/src/SURFEX/isba_to_topdsat.F90
+++ /dev/null
@@ -1,163 +0,0 @@
-!-------------------------------------------------------------------------------
-!     ####################
-      SUBROUTINE ISBA_TO_TOPDSAT(PKAPPA,PKAPPAC,KI,PRO_I,PRO_T)
-!     ####################
-!
-!!****  *ISBA_TO_TOPDSAT*  
-!!
-!!    PURPOSE
-!!    -------
-!     
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   23/11/2005
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_SURF_PAR,  ONLY : XUNDEF,NUNDEF
-!
-USE MODD_TOPODYN, ONLY : NNCAT, NNMC, NMESHT
-USE MODD_COUPLING_TOPD,ONLY: NMASKI, NMASKT, NNPIX
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN)              :: KI      ! Number of Isba meshes
-REAL, DIMENSION(:,:), INTENT(IN) :: PKAPPA  ! Hydrological indexes on the catchments 
-                                            ! at the previous time step
-REAL, DIMENSION(:), INTENT(IN)   :: PKAPPAC ! Hydrological index at saturation at the 
-                                            ! previous time step
-REAL, DIMENSION(:), INTENT(IN)   :: PRO_I   ! Runoff on Isba grid
-REAL, DIMENSION(:,:), INTENT(OUT):: PRO_T   ! Runoff on TOPODYN grid
-!
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                          :: JCAT, JPIX, JMESH_ISBA,JJ ! Loop indexes
-INTEGER, DIMENSION(KI)           :: INSAT       ! number of saturated pixels in an ISBA mesh
-INTEGER, DIMENSION(KI)           :: INDRY       ! Number of non-saturated pixels in an ISBA mesh
-REAL, DIMENSION(NNCAT,NMESHT)    :: ZROSAT      ! 
-REAL, DIMENSION(NNCAT,NMESHT)    :: ZRODRY      ! 
- CHARACTER(LEN=30)                :: YVAR        ! name of results file
-!
-REAL::ZSMALL,ZTMP,ZTMP2
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('ISBA_TO_TOPDSAT',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization :
-!               --------------
-!
-INSAT(:)=0
-INDRY(:)=0
-ZROSAT(:,:)=0.0
-ZRODRY(:,:)=0.0
-!
-! Only Isba meshes over studied catchments are scanned
-DO JMESH_ISBA = 1,KI
-  !
-  DO JCAT = 1,NNCAT
-    !
-    JJ=1
-    JPIX=NMASKI(JMESH_ISBA,JCAT,JJ)
-    !
-    DO WHILE (JPIX/=NUNDEF .AND.(JJ<=SIZE(NMASKI,3)))
-      !
-      IF (PKAPPA(JCAT,JPIX)/=XUNDEF .AND. NMASKT(JCAT,JPIX)/=NUNDEF) THEN
-        ! Calculation of the saturated and dry catchment pixels in each Isba mesh
-        IF (PKAPPA(JCAT,JPIX).GE.PKAPPAC(JCAT)) THEN
-          INSAT(NMASKT(JCAT,JPIX)) = INSAT(NMASKT(JCAT,JPIX)) + 1
-          ZROSAT(JCAT,JPIX) = PRO_I(NMASKT(JCAT,JPIX))
-        ELSE
-          INDRY(NMASKT(JCAT,JPIX)) = INDRY(NMASKT(JCAT,JPIX)) + 1 
-          ZRODRY(JCAT,JPIX) = PRO_I(NMASKT(JCAT,JPIX))
-        ENDIF
-      ENDIF
-      !
-      JJ=JJ+1
-      IF (JJ<=SIZE(NMASKI,3)) JPIX=NMASKI(JMESH_ISBA,JCAT,JJ)
-      !
-    ENDDO
-    !
-  ENDDO
-  !
-ENDDO
-!
-!
-DO JCAT = 1,NNCAT
-  !
-  DO JPIX = 1,NNMC(JCAT)
-    !
-    IF (NMASKT(JCAT,JPIX)/=NUNDEF) THEN
-      ! calculation of the runoff and deep drainage to rout in each Isba mesh, for each catchment
-      IF (INSAT(NMASKT(JCAT,JPIX)).GT.0 .AND. PKAPPA(JCAT,JPIX)/=XUNDEF) THEN
-        PRO_T(JCAT,JPIX) = ZROSAT(JCAT,JPIX) / INSAT(NMASKT(JCAT,JPIX))
-        ! if no runoff : calculation of the deep drainage to rout in each Isba mesh for each catchment
-      ELSEIF (INDRY(NMASKT(JCAT,JPIX)).GT.0 .AND. PKAPPA(JCAT,JPIX)/=XUNDEF) THEN
-        PRO_T(JCAT,JPIX) = ZRODRY(JCAT,JPIX) / INDRY(NMASKT(JCAT,JPIX))
-      ELSE
-        PRO_T(JCAT,JPIX) = 0.
-      ENDIF
-    ENDIF
-    !
-  ENDDO
-  !
-  ! budget control 
-  ZTMP=0.
-  ZTMP2=0.
-  !
-  DO JPIX = 1,NNMC(JCAT)
-    !
-    IF (PRO_T(JCAT,JPIX)/=XUNDEF) ZTMP = ZTMP + PRO_T(JCAT,JPIX)
-    IF ( NMASKT(JCAT,JPIX)/=NUNDEF) THEN
-      IF (PRO_I(NMASKT(JCAT,JPIX))/=XUNDEF .AND. NNPIX(NMASKT(JCAT,JPIX))/=0 ) &
-      ZTMP2 = ZTMP2 + PRO_I(NMASKT(JCAT,JPIX)) / NNPIX(NMASKT(JCAT,JPIX))
-    ENDIF
-    !
-  ENDDO!JPIX
-  !
-  ZSMALL=ABS(ZTMP2*0.001)
-  !
-  IF( ABS(ZTMP-ZTMP2) > ZSMALL ) THEN
-    WHERE ( PRO_T(JCAT,:)/=XUNDEF )
-      PRO_T(JCAT,:) = PRO_T(JCAT,:)- ((ZTMP-ZTMP2)/NNMC(JCAT))
-    ENDWHERE
-  ENDIF
-  !
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('ISBA_TO_TOPDSAT',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE ISBA_TO_TOPDSAT
diff --git a/src/SURFEX/make_mask_isba_to_topd.F90 b/src/SURFEX/make_mask_isba_to_topd.F90
deleted file mode 100644
index ced694a9facdd237b9b5c3b69f24ec0175bfb040..0000000000000000000000000000000000000000
--- a/src/SURFEX/make_mask_isba_to_topd.F90
+++ /dev/null
@@ -1,104 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE MAKE_MASK_ISBA_TO_TOPD(KI)
-!     #######################
-!
-!!****  *MAKE_MASK_ISBA_TO_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!
-!     Create a mask for each Surfex mesh and each catchment. 
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!                          
-!!    REFERENCE
-!!    ---------
-     
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* CNRM * 
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   16/03/2005
-!!                 11/2011 : Loops simplified (Vincendon)
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, NNMC
-USE MODD_COUPLING_TOPD, ONLY : NMASKT, NMASKI, NNPIX
-USE MODD_SURF_PAR,        ONLY : NUNDEF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN) :: KI    ! Grid dimensions
-!
-!*      0.2    declarations of local variables
-!
-INTEGER, DIMENSION(KI)  :: INBPIX_IN_MESH ! number of pixel in each ISBA mesh
-INTEGER                 :: JCAT, JPIX, INUMPIX
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_ISBA_TO_TOPD',0,ZHOOK_HANDLE)
-!
-INUMPIX=MAXVAL(NNPIX)
-!
-ALLOCATE(NMASKI(KI,NNCAT,INUMPIX))
-NMASKI(:,:,:) = NUNDEF
-!
-INBPIX_IN_MESH(:) = 0
-!
-DO JCAT = 1,NNCAT
-  !
-  DO JPIX = 1,NNMC(JCAT)
-    !si le point du bassin versant est dans une maille isba
-    IF ((NMASKT(JCAT,JPIX)/=0).AND.(NMASKT(JCAT,JPIX)/=NUNDEF)) THEN
-      !indice du point du bassin versant dans la maille isba
-      INBPIX_IN_MESH(NMASKT(JCAT,JPIX)) = INBPIX_IN_MESH(NMASKT(JCAT,JPIX)) + 1
-      ! nmaski associe à la maille isba, au bassin versant et au numÊro du point 
-      ! du bassin versant dans la maille isba, l'indice du point dans le bassin
-      ! versant
-      NMASKI(NMASKT(JCAT,JPIX),JCAT,INBPIX_IN_MESH(NMASKT(JCAT,JPIX))) = JPIX
-    ENDIF
-    !
-  ENDDO
-  !
-ENDDO
-! write(*,*) 'NMASKT min et max',MINVAL(NMASKT(1,:)),MAXVAL(NMASKT(1,:))
-! write(*,*) 'NMASKT min et max',MINVAL(NMASKT(2,:)),MAXVAL(NMASKT(2,:))
-! write(*,*) 'NMASKT min et max',MINVAL(NMASKT(3,:)),MAXVAL(NMASKT(3,:))
-! write(*,*) 'NMASKT min et max',MINVAL(NMASKT(4,:)),MAXVAL(NMASKT(4,:))
-! write(*,*) 'NMASKI 3132 min et max',MINVAL(NMASKI(MINVAL(NMASKT(1,:)),1,:)),MAXVAL(NMASKI(MINVAL(NMASKT(1,:)),1,:))
-! write(*,*) 'NMASKI 6662 min et max',MINVAL(NMASKI(MAXVAL(NMASKT(1,:)),1,:)),MAXVAL(NMASKI(MAXVAL(NMASKT(1,:)),1,:))
-! stop
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_ISBA_TO_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE MAKE_MASK_ISBA_TO_TOPD
-
-
-
-
-
-
-
diff --git a/src/SURFEX/make_mask_topd_to_isba.F90 b/src/SURFEX/make_mask_topd_to_isba.F90
deleted file mode 100644
index e3b117d0d56c1cbdabb8c82d155da556ebaa4901..0000000000000000000000000000000000000000
--- a/src/SURFEX/make_mask_topd_to_isba.F90
+++ /dev/null
@@ -1,230 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-    SUBROUTINE MAKE_MASK_TOPD_TO_ISBA(KI)
-!     #######################
-!
-!!****  *MAKE_MASK_TOPD_TO_ISBA(*  
-!!
-!!    PURPOSE
-!!    -------
-!
-!     Create a mask for each catchment. 
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!                     
-!!    REFERENCE
-!!    ---------
-     
-!!    AUTHOR
-!!    ------
-!!
-!!      L. Labatut & K. Chancibault	* CNRM *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   16/03/2005
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN,       ONLY: NNCAT, NNYC, XY0, XDXT, NNXC,&
-                                XX0, XTOPD, XNUL, NLINE, NMESHT
-USE MODD_COUPLING_TOPD, ONLY: NMASKT
-USE MODD_SURF_PAR,        ONLY : XUNDEF, NUNDEF
-!
-USE MODI_GET_LUOUT
-USE MODI_ABOR1_SFX
-USE MODI_WRITE_FILE_MAP
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN) :: KI    ! Grid dimensions
-!
-!*      0.2    declarations of local variables
-!
- CHARACTER(LEN=30)  :: YVAR        ! name of results file
-INTEGER            :: JCAT, JJ, JI, IDX     ! loop control  
-INTEGER            :: II,ILINE ! work integer variables
-INTEGER            :: IDXM     ! indexes of Isba grid meshes and nodes
-INTEGER            :: ILUOUT   ! unit
-REAL               :: ZXT, ZYT ! catchment grid nodes Lambert II coordinates 
-REAL               :: ZX1, ZX2, ZX3, ZX4, ZY1, ZY2, ZY3, ZY4  ! Isba mesh Lambert II coordinates
-REAL               :: ZXA, ZXB, ZYA, ZYB
-REAL, DIMENSION(NNCAT,NMESHT):: ZWRK
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT('OFFLIN',ILUOUT)
-!  
-!WRITE(*,*) 'Y0=', XY0(1), XY0(2)
-!WRITE(*,*) 'X0=', XX0(1), XX0(2)
-!WRITE(*,*) 'X et Y isba=', XXI(1), XXI(180), XYI(1), XYI(180)
-!loop on catchments
-DO JCAT=1,NNCAT
-  !
-  IDXM  = 1
-  CALL INIT_4POINTS(IDXM,ZX1,ZX2,ZX3,ZX4,ZY1,ZY2,ZY3,ZY4)
-  !
-  DO JJ=1,NNYC(JCAT)   ! number of topographic points on the y axis
-    !
-    ZYT = XY0(JCAT) + (JJ-1) * XDXT(JCAT)
-    ZYT = ZYT + 0.5 * XDXT(JCAT)
-    !
-    DO JI=1,NNXC(JCAT)
-      !
-      ZXT = XX0(JCAT) + (JI-1) * XDXT(JCAT)
-      ZXT = ZXT + 0.5 * XDXT(JCAT)
-      !
-      IDX = (JJ-1) * NNXC(JCAT) + JI   ! index of the point among all in the catchment
-      !
-      !* on vérifie que le pixel du MNT appartient au BV
-      IF ( XTOPD(JCAT,IDX).NE.XNUL(JCAT) ) THEN
-        !* calcule des coordonées X et Y inf et sup de la maille ISBA considérée
-        CALL GET_COORD(ZXT,ZYT,ZX1,ZX2,ZX3,ZX4,ZY1,ZY2,ZY3,ZY4,ZXA,ZYA,ZXB,ZYB)
-        !* si on se trouve sur le premier pixel du MNT ou si le pixel du MNT n'est pas 
-        !dans la maille Isba considérée (qui est celle dans laquelle se trouve le pixel précédent)
-        IF (ZXT.LT.ZXA.OR.ZXT.GE.ZXB.OR.ZYT.LT.ZYA.OR.ZYT.GE.ZYB) THEN
-          !* on repart de la premičre maille de la grille Isba
-          IDXM = 1
-          CALL INIT_4POINTS(IDXM,ZX1,ZX2,ZX3,ZX4,ZY1,ZY2,ZY3,ZY4)
-          CALL GET_COORD(ZXT,ZYT,ZX1,ZX2,ZX3,ZX4,ZY1,ZY2,ZY3,ZY4,ZXA,ZYA,ZXB,ZYB)
-          !* on parcours les mailles de la grille Isba, jusqu'ā ce qu'on trouve la maille ā laquelle appartient le pixel du MNT
-          DO WHILE (ZXT.LT.ZXA.OR.ZXT.GE.ZXB.OR.ZYT.LT.ZYA.OR.ZYT.GE.ZYB)
-            IDXM = IDXM + 1
-            IF (IDXM.GE.KI) THEN
-              WRITE(*,*) 'ZXT', ZXT,'ZYT',ZYT
-              WRITE(*,*) 'indices Isba:',IDXM,'>=',KI
-              CALL ABOR1_SFX("MAKE_MASK_TOPD_TO_ISBA: PROBLEM")
-            ENDIF
-            CALL INIT_4POINTS(IDXM,ZX1,ZX2,ZX3,ZX4,ZY1,ZY2,ZY3,ZY4)
-            CALL GET_COORD(ZXT,ZYT,ZX1,ZX2,ZX3,ZX4,ZY1,ZY2,ZY3,ZY4,ZXA,ZYA,ZXB,ZYB)
-          ENDDO
-        ENDIF
-        IF (NLINE(JCAT,IDX)/=0) NMASKT(JCAT,NLINE(JCAT,IDX)) = IDXM
-      ENDIF
-    ENDDO
-  ENDDO
-ENDDO
-!
-YVAR='.mask_topd'
-WHERE (NMASKT(:,:)/=NUNDEF)
-  ZWRK(:,:)=REAL(NMASKT)
-ELSEWHERE
-  ZWRK(:,:)=XUNDEF
-ENDWHERE
- CALL WRITE_FILE_MAP(ZWRK,YVAR)
-!CALL WRITE_FILE_MAP(REAL(NMASKT),YVAR)
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA',1,ZHOOK_HANDLE)
-!
-CONTAINS
-!
-SUBROUTINE INIT_4POINTS(KDXM,PX1,PX2,PX3,PX4,PY1,PY2,PY3,PY4)
-!
-USE MODD_COUPLING_TOPD, ONLY: NIMAX, XXI, XYI
-!
-INTEGER, INTENT(IN) :: KDXM
-REAL, INTENT(OUT) :: PX1, PX2, PX3, PX4
-REAL, INTENT(OUT) :: PY1, PY2, PY3, PY4
-!
-INTEGER :: ILINE, II, IDXN
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA:INIT_4POINTS',0,ZHOOK_HANDLE)
-!
-ILINE = INT(KDXM/(NIMAX))+1      ! number of the current line
-II    = KDXM-((ILINE-1)*NIMAX)   ! index of point in the line
-IDXN  = (ILINE-1)*(NIMAX+1)+II   ! indice du point dans la grille isba
-!
-PX1 = XXI(IDXN)              ! coordonnée X du point courant
-PX2 = XXI(IDXN+1)            ! coordonnée X du point suivant
-PX3 = XXI(IDXN+(NIMAX+1))    ! coordonnée X du point aligné sur la ligne suivante 
-PX4 = XXI(IDXN+1+(NIMAX+1))  ! coordonnée X du point suivant sur la ligne suivante
-!
-PY1 = XYI(IDXN)              ! coordonnée Y du point courant
-PY2 = XYI(IDXN+1)            ! coordonnée Y du point suivant
-PY3 = XYI(IDXN+(NIMAX+1))    ! coordonnée Y du point aligné sur la ligne suivante
-PY4 = XYI(IDXN+1+(NIMAX+1))  ! coordonnée Y du point suivant sur la ligne suivante
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA:INIT_4POINTS',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE INIT_4POINTS
-!
-SUBROUTINE GET_COORD(PXT,PYT,PX1,PX2,PX3,PX4,PY1,PY2,PY3,PY4,&
-                     PXA,PYA,PXB,PYB)
-!
-REAL, INTENT(IN) :: PXT, PYT
-REAL, INTENT(IN) :: PX1, PX2, PX3, PX4
-REAL, INTENT(IN) :: PY1, PY2, PY3, PY4
-REAL, INTENT(OUT) :: PXA, PYA, PXB, PYB
-REAL :: ZFA, ZFB
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA:GET_COORD',0,ZHOOK_HANDLE)
-!
-IF((PX3-PX1).EQ.0.0) THEN
-  PXA = PX1
-ELSE
-  CALL GET_LINE_PARAM(PX1,PY1,PX3,PY3,ZFA,ZFB)
-  PXA = (ZYT-ZFB)/ZFA
-ENDIF
-!
-IF ((PX4-PX2).EQ.0.0) THEN
-  PXB = PX2
-ELSE
-  CALL GET_LINE_PARAM(PX2,PY2,PX4,PY4,ZFA,ZFB)
-  PXB = (ZYT-ZFB)/ZFA
-ENDIF
-!
-IF ((PY2-PY1).EQ.0.0) THEN
-  PYA = PY2
-ELSE
-  CALL GET_LINE_PARAM(PX1,PY1,PX2,PY2,ZFA,ZFB)
-  PYA = ZFA*ZXT+ZFB
-ENDIF
-!
-IF ((PY4-PY3).EQ.0.0) THEN
-  PYB = PY4
-ELSE
-  CALL GET_LINE_PARAM(PX3,PY3,PX4,PY4,ZFA,ZFB)
-  PYB = ZFA*ZXT+ZFB
-ENDIF
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA:GET_COORD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE GET_COORD
-!
-SUBROUTINE GET_LINE_PARAM(PX1,PY1,PX2,PY2,PFA,PFB)
-!
-REAL, INTENT(IN) :: PX1, PX2, PY1, PY2
-REAL, INTENT (OUT) :: PFA, PFB
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA:GET_LINE_PARAM',0,ZHOOK_HANDLE)
-!
-PFA = (PY2 - PY1) / (PX2 - PX1)  ! slope
-PFB = PY1 - PFA * PX1 ! offset
-!
-IF (LHOOK) CALL DR_HOOK('MAKE_MASK_TOPD_TO_ISBA:GET_LINE_PARAM',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE GET_LINE_PARAM
-!
-END SUBROUTINE MAKE_MASK_TOPD_TO_ISBA
diff --git a/src/SURFEX/modd_budget_coupl_rout.F90 b/src/SURFEX/modd_budget_coupl_rout.F90
deleted file mode 100644
index 9d5b7f1d109b5fea04f8cc56033a0c4f3733593c..0000000000000000000000000000000000000000
--- a/src/SURFEX/modd_budget_coupl_rout.F90
+++ /dev/null
@@ -1,81 +0,0 @@
-!     ######spl
-      MODULE MODD_BUDGET_COUPL_ROUT
-!     ######################
-!
-!!****  *MODD_BUDGET_COUPL_ROUT* - declaration of variables
-!                                  useful for budget computations when
-!                                  coupling with TOPMODEL
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!**  IMPLICIT ARGUMENTS
-!!    ------------------
-!!      None 
-!!
-!!    REFERENCE
-!!    ---------
-!!                
-!!    AUTHOR
-!!    ------
-!!	B. Vincendon *Meteo France*
-!!
-!!    MODIFICATIONS
-!!    -------------
-!-------------------------------------------------------------------------------
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-IMPLICIT NONE
-!
-! Water entering the system
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_RAIN, XB_SNOW ! Rain and Snow
-! Water going out of the system
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_WR ! Interception
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_EVAP! Evaporation
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_RUNOFF_TOPD,XB_RUNOFF_ISBA !Runoff
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_DRAIN !Drainage
-! Water in the ground
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_WG2, XB_WG3, XB_WGTOT !liquid
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_WGI2, XB_WGI3, XB_WGITOT !solid
-!REAL,ALLOCATABLE,DIMENSION(:) :: XB_DWR_TOT,XB_DWI_TOT, XB_DW_TOT
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_SWE1, XB_SWE2, XB_SWE3, XB_SWETOT ! snow melt
-!REAL,ALLOCATABLE,DIMENSION(:) :: XB_RUN2,XB_RUN3!bv pour verif
-!
-! Values at the previous time step
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_WRM
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_EVAPM
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_RUNOFF_TOPDM,XB_RUNOFF_ISBAM 
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_DRAINM
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_WG2M, XB_WG3M, XB_WGTOTM
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_WGI2M,XB_WGI3M, XB_WGITOTM
-REAL,ALLOCATABLE,DIMENSION(:) :: XB_SWE1M, XB_SWE2M, XB_SWE3M, XB_SWETOTM
-!REAL,ALLOCATABLE,DIMENSION(:) :: XB_RUN2M,XB_RUN3M!bv pour verif
-! Useful
-REAL,ALLOCATABLE,DIMENSION(:)   :: XB_DG2, XB_DG3
-REAL,ALLOCATABLE,DIMENSION(:)   :: XB_MESH_SIZE
-REAL,ALLOCATABLE,DIMENSION(:,:) :: XB_ABV_BYMESH !fraction de chaque BV ds chaque maille 
-REAL,ALLOCATABLE,DIMENSION(:)   :: XB_ATOP_BYMESH!fraction de tous BV ds chaque maille 
-!
-! Variable to keep and write budget terms
-REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XB_VAR_BV
-REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XB_VAR_NOBV
-REAL,ALLOCATABLE,  DIMENSION(:,:)   :: XB_VAR_TOT
-REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XB_VAR_Q
-!
- CHARACTER(LEN=6),ALLOCATABLE,  DIMENSION(:)     :: YB_VAR  
-!
-! Discharge variables
-!
-REAL, ALLOCATABLE, DIMENSION(:)  :: XB_QTOT,XB_QRUN
-REAL, ALLOCATABLE, DIMENSION(:)  :: XB_QDR
-REAL, ALLOCATABLE, DIMENSION(:)  :: XB_STOCK_TOT,XB_STOCK_RUN
-REAL, ALLOCATABLE, DIMENSION(:)  :: XB_STOCK_DR
-! Values at the previous time step
-REAL, ALLOCATABLE, DIMENSION(:)  :: XB_QTOTM,XB_QRUNM
-REAL, ALLOCATABLE, DIMENSION(:)  :: XB_QDRM
-!
- CHARACTER(LEN=6),ALLOCATABLE,  DIMENSION(:)  :: YB_VARQ  
-!
-END MODULE MODD_BUDGET_COUPL_ROUT
diff --git a/src/SURFEX/modd_coupling_topd.F90 b/src/SURFEX/modd_coupling_topd.F90
deleted file mode 100644
index 89b833dccc989755f2fd91f1942b80d9a6c45266..0000000000000000000000000000000000000000
--- a/src/SURFEX/modd_coupling_topd.F90
+++ /dev/null
@@ -1,116 +0,0 @@
-!     ###########################
-      MODULE MODD_COUPLING_TOPD
-!     ###########################
-!
-!!****  *MODD_COUPLING_TOPD - declaration of exchanged variables from Topodyn to ISBA
-!!
-!!    PURPOSE
-!!    -------
-!
-!!
-!!**  IMPLICIT ARGUMENTS
-!!    ------------------
-!!      None 
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!     F. Habets and K. Chancibault
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original       29/09/03
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-IMPLICIT NONE
-!
-!-------------------------------------------------------------------------------
-!
-LOGICAL                           :: LCOUPL_TOPD      !if T, performs coupling with Topmodel
-LOGICAL                           :: LBUDGET_TOPD     !if T, computes budget
-LOGICAL                           :: LTOPD_STEP
-!
-INTEGER                             :: NTOPD_STEP
-INTEGER                             :: NFREQ_MAPS_WG       !frequency of output WG maps
-INTEGER                             :: NFREQ_MAPS_ASAT     !frequency of output ASAT maps
-!
-INTEGER                             :: NNB_TOPD   ! Ratio between Time steps of Topmodel and ISBA
-!
-INTEGER                             :: NIMAX     ! number of ISBA grid points on 
-                                                 ! abscissa axis
-INTEGER                             :: NJMAX     ! number of ISBA grid points on ordinate 
-                                                 ! axis
-REAL, ALLOCATABLE, DIMENSION(:)     :: XXI ! Extended Lambert II coordinates of Isba 
-REAL, ALLOCATABLE, DIMENSION(:)     :: XYI ! nodes 
-!
-INTEGER, ALLOCATABLE, DIMENSION(:)  :: NNPIX     ! Number of Topmodel pixels in an ISBA mesh
-INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: NMASKI ! pixel number of each catchment in each isba mesh
-INTEGER, ALLOCATABLE, DIMENSION(:,:):: NMASKT    ! mask
-!
-REAL, ALLOCATABLE, DIMENSION(:)     :: XAS_NATURE ! Packed contributive area fraction on Nature grid
-REAL, ALLOCATABLE, DIMENSION(:)     :: XATOP      ! Packed area fraction WITH TOPMODEL on Nature grid
-!
-INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: NNBV_IN_MESH   ! Number of pixel of a partical cathment in an ISBA mesh
-REAL, ALLOCATABLE, DIMENSION(:,:)     :: XBV_IN_MESH    ! Area of the ISBA meshes covered by a partical cathment
-REAL, ALLOCATABLE, DIMENSION(:)       :: XTOTBV_IN_MESH ! Area of the ISBA meshes covered by all cathments
-!
-REAL, ALLOCATABLE, DIMENSION(:)     :: XDTOPI   ! depth of the soil for lateral 
-                                                ! distribution on ISBA grid (m)
-REAL, ALLOCATABLE, DIMENSION(:,:)   :: XDTOPT   ! depth of the Isba soil on TOP-LAT 
-                                                ! grid (m)
-!
-REAL, ALLOCATABLE, DIMENSION(:)     :: XWG_FULL   ! Water content from Isba on the full domain
-REAL, ALLOCATABLE, DIMENSION(:,:)   :: XWGT     ! ISBA water content 
-!
-REAL, ALLOCATABLE, DIMENSION(:)     :: XWSTOPI  ! total water content at saturation (m3/m3)
-                                                ! on XDTOPI on ISBA grid
-REAL, ALLOCATABLE, DIMENSION(:,:)   :: XWSTOPT  ! total water content at saturation (m3/m3)
-                                                ! on XDTOPT on TOP-LAT grid
-REAL, ALLOCATABLE, DIMENSION(:)     :: XWFCTOPI ! total field capacity on XDTOPI (m3/m3)
-REAL, ALLOCATABLE, DIMENSION(:,:)   :: XWFCTOPT ! total field capacity on XDTOPT (m3/m3)
-REAL, ALLOCATABLE, DIMENSION(:)     :: XCSTOPI  ! hydraulic conductivity at saturation on 
-                                                ! Isba grid, on XDTOPI
-REAL, ALLOCATABLE, DIMENSION(:,:)   :: XWTOPT   ! water storage on TOP-LAT grid, after
-                                                ! lateral distribution
-!
-! * pour bilans
-REAL, ALLOCATABLE, DIMENSION(:)       :: XAVG_RUNOFFCM !cumulated runoff  (kg/m2) at t-dt
-REAL, ALLOCATABLE, DIMENSION(:)       :: XAVG_DRAINCM ! cumulated drainage calculated from Isba (kg/m2) at t-dt
-!
-REAL, ALLOCATABLE, DIMENSION(:,:)     :: XKA_PRE   ! Hydrological indexes at the previous time step
-REAL, ALLOCATABLE, DIMENSION(:)       :: XKAC_PRE  ! Hydrological index at saturation at the previous time step
-!
-REAL, ALLOCATABLE, DIMENSION(:,:)     :: XDMAXFC   ! Deficit at the field capacity level
-REAL, ALLOCATABLE, DIMENSION(:)       :: XWSUPSAT  ! pour calculer le volume d'eau perdu au-dessus de la saturation
-!
-REAL, ALLOCATABLE, DIMENSION(:)      ::  XDRAIN_TOP ! Value of drainage on TOPMODEL grid
-REAL, ALLOCATABLE, DIMENSION(:)      ::  XRUNOFF_TOP! Value of runoff on TOPMODEL grid
-!
-REAL, ALLOCATABLE, DIMENSION(:)      ::  XFRAC_D2 ! fraction of the second layer concerned with lateral transferts
-REAL, ALLOCATABLE, DIMENSION(:)      ::  XFRAC_D3 ! fraction of the third layer concerned with lateral transferts
-!
-REAL, ALLOCATABLE, DIMENSION(:)      ::  XWGI_FULL ! soil ice content
-!
-REAL, ALLOCATABLE, DIMENSION(:,:)    :: XRUN_TOROUT,XDR_TOROUT
-!
-LOGICAL                              :: LSTOCK_TOPD ! true to stock runoff and drainage values (for another simulation)
-!
-INTEGER                              :: NNB_STP_RESTART ! number of time step to restart from a previous simulation
-INTEGER                              :: NNB_STP_STOCK   ! number of time step to write for the next simulation
-!
-INTEGER, DIMENSION(:), ALLOCATABLE :: NYEAR      ! Year of the beginning of the simulation.
-INTEGER, DIMENSION(:), ALLOCATABLE :: NMONTH     ! Month of the beginning of the simulation.
-INTEGER, DIMENSION(:), ALLOCATABLE :: NDAY      ! Date of the beginning of the simulation.
-INTEGER, DIMENSION(:), ALLOCATABLE :: NH      ! Hour of the beginning of the simulation.
-INTEGER, DIMENSION(:), ALLOCATABLE :: NM      ! Minutes of the beginning of the simulation.
-!
-! **** For special f, dc exponential profile
-REAL, DIMENSION(:), ALLOCATABLE :: XF_PARAM
-REAL, DIMENSION(:), ALLOCATABLE :: XC_DEPTH_RATIO
-!
-END MODULE MODD_COUPLING_TOPD
-
diff --git a/src/SURFEX/modd_dummy_exp_profile.F90 b/src/SURFEX/modd_dummy_exp_profile.F90
deleted file mode 100644
index e071623fbb8e68f7c20f5fb9b26f6e0fb7c35f0b..0000000000000000000000000000000000000000
--- a/src/SURFEX/modd_dummy_exp_profile.F90
+++ /dev/null
@@ -1,43 +0,0 @@
-!     ###########################
-      MODULE MODD_DUMMY_EXP_PROFILE
-!     ###########################
-!
-!!****  *MODD_DUMMY_EXP_PROFILE - declaration For special f, dc exponential profile
-!!
-!!    PURPOSE
-!!    -------
-!
-!!
-!!**  IMPLICIT ARGUMENTS
-!!    ------------------
-!!      None 
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!     B. Vincendon
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original       10/08/11
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-USE MODD_TOPD_PAR, ONLY : JPCAT
-!
-IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! **** For special f, dc exponential profile
-! values for each isba_mesh
-REAL, DIMENSION(:), ALLOCATABLE :: XF_PARAM
-REAL, DIMENSION(:), ALLOCATABLE :: XC_DEPTH_RATIO
-!values for each catchment
-REAL, DIMENSION(JPCAT) :: XF_PARAM_BV
-REAL, DIMENSION(JPCAT) :: XC_DEPTH_RATIO_BV
-!-------------------------------------------------------------------------------------
-!
-END MODULE MODD_DUMMY_EXP_PROFILE
-
diff --git a/src/SURFEX/modd_topd_par.F90 b/src/SURFEX/modd_topd_par.F90
deleted file mode 100644
index 4b45f18385c4e0eb010e416ffd5983fa2fae7ad9..0000000000000000000000000000000000000000
--- a/src/SURFEX/modd_topd_par.F90
+++ /dev/null
@@ -1,48 +0,0 @@
-!-------------------------------------------------------------------------------
-!     ######################
-      MODULE MODD_TOPD_PAR
-!     ######################
-!
-!!****  *MODD_TOPD_PAR* - declaration of parameter variables
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this declarative module is to specify  the variables 
-!     which have the PARAMETER attribute   
-!
-!!
-!!**  IMPLICIT ARGUMENTS
-!!    ------------------
-!!      None 
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book2 of documentation of Meso-NH (module MODD_PARAMETER)
-!!          
-!!    AUTHOR
-!!    ------
-!!	V. Ducrocq   *Meteo France*
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    4/07/94                      
-!!      Modification 10/03/95 (I.Mallet)   add the coupling files maximum number
-!!      Modification 10/04/95 (Ph. Hereil) add the budget related informations
-!-------------------------------------------------------------------------------
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-IMPLICIT NONE
-!
-REAL,    PARAMETER :: XSTEPK = 0.05   ! discretization step of the saturation 
-                                      ! index KAPPA
-INTEGER, PARAMETER :: NDIM = 20       ! dimension of the XCONN array third 
-                                      ! index 
-INTEGER, PARAMETER :: JPCAT = 10      ! number max of catchments
-!
-!values for each catchment
-REAL, DIMENSION(JPCAT) :: XF_PARAM_BV
-REAL, DIMENSION(JPCAT) :: XC_DEPTH_RATIO_BV
-!
-END MODULE MODD_TOPD_PAR
diff --git a/src/SURFEX/pgd_topd.F90 b/src/SURFEX/pgd_topd.F90
deleted file mode 100644
index 4588544ee388af6e3ff2ce850f531eb090b7cce4..0000000000000000000000000000000000000000
--- a/src/SURFEX/pgd_topd.F90
+++ /dev/null
@@ -1,360 +0,0 @@
-!-------------------------------------------------------------------------------
-!     #############################################################
-      SUBROUTINE PGD_TOPD(HPROGRAM)
-!     #############################################################
-!
-!!****  *PGD_TOPD* - routine to determine the masks that permit to couple ISBA grid with Topmodel one
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!**  METHOD
-!!    ------
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!    REFERENCE
-!!    ---------
-!!   (from initial version of init_coupl_topo.f90 by K. Chancibault, modified by
-!!     M. Lelay and B. Vincendon)
-!!
-!!    AUTHOR
-!!    ------
-!!	B. Vincendon   *Meteo France*	
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    11/2011
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-USE MODD_TOPODYN,         ONLY : CCAT, NNCAT, XRTOP_D2, NMESHT, XDXT
-USE MODD_COUPLING_TOPD,   ONLY : LCOUPL_TOPD, NIMAX, NJMAX, &
-                                 XXI, XYI, NMASKI, NMASKT, NNPIX,&
-                                 NNBV_IN_MESH, XBV_IN_MESH, XTOTBV_IN_MESH
-USE MODD_DUMMY_EXP_PROFILE, ONLY : XF_PARAM_BV, XC_DEPTH_RATIO_BV
-!
-USE MODD_SURF_PAR,          ONLY : NUNDEF
-USE MODD_SURF_ATM_GRID_N,   ONLY : XGRID_PAR, CGRID  !mll: ajout CGRID
-USE MODD_SURF_ATM_n,        ONLY : NDIM_FULL
-!
-USE MODD_ISBA_n, ONLY : CISBA
-!
-USE MODE_GRIDTYPE_CONF_PROJ
-USE MODE_GRIDTYPE_LONLAT_REG
-USE MODE_GRIDTYPE_IGN
-!
-USE MODI_GET_LUOUT
-USE MODI_READ_NAM_PGD_TOPD
-USE MODI_INIT_TOPD
-USE MODI_ABOR1_SFX
-USE MODI_MAKE_MASK_TOPD_TO_ISBA
-USE MODI_MAKE_MASK_ISBA_TO_TOPD
-USE MODI_WRITE_FILE_MASKTOPD
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-USE MODI_TOPD_TO_ISBA_SLOPE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*       0.1   Declarations of arguments
-!              -------------------------
-!
- CHARACTER(LEN=*),  INTENT(IN)     :: HPROGRAM    !
-!
- CHARACTER(LEN=50),DIMENSION(NNCAT) :: CNAME
-INTEGER                   :: IL                     ! number of points
-INTEGER                   :: JJ,JI,JK,JWRK ! loop control 
-INTEGER                   :: JCAT,JMESH,JPIX ! loop control 
-INTEGER                           :: ILUOUT       ! Logical unit for output filr
-INTEGER                           :: IUNIT       ! Logical unit for carte_asat file
-INTEGER                           :: IMESHL       !  number of ISBA grid nodes
-!
-REAL, DIMENSION(:), ALLOCATABLE   :: ZXI, ZYI     ! natural coordinates of ISBA grid (conformal projection or latlon)
-REAL, DIMENSION(:), ALLOCATABLE   :: ZDXI, ZDYI   ! Isba grid resolution in the conformal projection
-REAL, DIMENSION(:), ALLOCATABLE   :: ZXN, ZYN     ! isba nodes coordinates in the Lambert II coordinates
-REAL, DIMENSION(:), ALLOCATABLE   :: ZLAT,ZLON    ! Isba nodes geographical coordinates
-REAL, DIMENSION(:), ALLOCATABLE   :: ZDTAV        ! Averaged depth soil on TOP-LAT grid
-REAL                              :: ZLAT0    ! reference latitude
-REAL                              :: ZLON0    ! reference longitude
-REAL                              :: ZLONMIN,ZLONMAX  ! min and max longitude values (latlon coordinates)
-REAL                              :: ZLATMIN,ZLATMAX  ! min and max latitude values (latlon coordinates)
-REAL                              :: ZRPK     ! projection parameter 
-!                                             !   K=1 : stereographic north pole
-!                                             ! 0<K<1 : Lambert, north hemisphere
-!                                             !   K=0 : Mercator
-!                                             !-1<K<0 : Lambert, south hemisphere
-!                                             !   K=-1: stereographic south pole
-REAL                              :: ZBETA    ! angle between grid and reference longitude
-REAL                              :: ZLATOR   ! latitude  of point of coordinates X=0, Y=0
-REAL                              :: ZLONOR   ! longitude of point of coordinates X=0, Y=0!
-REAL, DIMENSION(:), ALLOCATABLE   :: ZF_PARAM,ZC_DEPTH_RATIO  !
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('PGD_TOPD',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!  
- CALL READ_NAM_PGD_TOPD(HPROGRAM,LCOUPL_TOPD,CCAT,XF_PARAM_BV,XC_DEPTH_RATIO_BV)
-!
-IF (LCOUPL_TOPD .AND. CISBA/='3-L') &
-  CALL ABOR1_SFX("PGD_TOPD: coupling with topmodel only runs with CISBA=3-L")
-!
-!         1.   Reads the namelists
-!              --------------------
-IF (LCOUPL_TOPD) THEN
-  !
-  WRITE(ILUOUT,*) 'Debut pgd_topd'
-  !
-  !         3.   Initialises variables specific to Topmodel
-  !              -------------------------------------------
-  WRITE(ILUOUT,*) 'NNCAT',NNCAT
-  !
-  CALL INIT_TOPD(HPROGRAM)
-  !
-  !         4.   Compute masks to couple ISBA and TOPODYN grids
-  !              -------------------------------------------------------
-  !
-  !*      1.      Calcul of the plane coordinates 
-  !               -------------------------------
-  !
-  !*      1.1 Gestion des projections conformes
-  !           ---------------------------------
-  !
-  ALLOCATE(NMASKT(NNCAT,NMESHT))
-  NMASKT(:,:)=NUNDEF
-  !
-  IF(CGRID.EQ.'CONF PROJ') THEN
-    !
-    WRITE(ILUOUT,*) 'GRILLE PROJ CONF (application Cevennes)'
-    !
-    !*      1.1.1   lecture des coordonnees X et Y conformes 
-    !               -----------------------------------------
-    !
-    ALLOCATE(ZXI(NDIM_FULL))
-    ALLOCATE(ZYI(NDIM_FULL))
-    ZXI(:)=0.0
-    ZYI(:)=0.0
-    !
-    ALLOCATE(ZDXI(NDIM_FULL))
-    ALLOCATE(ZDYI(NDIM_FULL))
-    !
-    CALL GET_GRIDTYPE_CONF_PROJ(XGRID_PAR,PLAT0=ZLAT0,PLON0=ZLON0,PRPK=ZRPK, &
-                                PBETA=ZBETA,PLATOR=ZLATOR,PLONOR=ZLONOR,     &
-                                KIMAX=NIMAX,KJMAX=NJMAX,PX=ZXI,PY=ZYI,       &
-                                PDX=ZDXI,PDY=ZDYI)
-    !
-    IMESHL = (NIMAX+1)*(NJMAX+1)
-    !
-    !*      1.1.2   calcul des coordonnees X et Y conformes des noeuds de la grille ISBA
-    !               --------------------------------------------------------------------
-    ALLOCATE(ZXN(IMESHL))
-    ALLOCATE(ZYN(IMESHL))
-    !
-    DO JJ=1,NIMAX
-      ZXN(JJ) = ZXI(JJ) - ZDXI(JJ)/2.
-    ENDDO
-    ZXN(NIMAX+1) = ZXI(NIMAX) + ZDXI(NIMAX)/2.
-    !
-    DO JJ=1,NJMAX
-      JWRK = (JJ-1)*(NIMAX+1)+1   ! indice sur la grille des noeuds
-      JI = (JJ-1)*NIMAX+1         ! indice sur la grille des mailles
-      ZYN(JWRK) = ZYI(JI) - ZDYI(JI)/2.
-    ENDDO
-    !
-    JJ = ((NJMAX+1)-1)*(NIMAX+1)+1  ! Indice sur la grille des noeuds
-    JI = (NJMAX-1)*NIMAX+1
-    ZYN(JJ) = ZYI(JI) + ZDYI(JI)/2.
-    !
-    DEALLOCATE(ZDXI)
-    DEALLOCATE(ZDYI)
-    DEALLOCATE(ZXI)
-    DEALLOCATE(ZYI)
-    !
-    DO JJ=1,NIMAX+1
-      DO JI=2,NJMAX+1
-        JK = (JI-1)*(NIMAX+1)+JJ
-        ZXN(JK) = ZXN(JJ)
-      ENDDO
-    ENDDO
-    !
-    DO JI=1,NJMAX+1
-      DO JJ=2,NIMAX+1
-        JK = (JI-1)*(NIMAX+1)+JJ
-        JWRK = (JI-1)*(NIMAX+1)+1
-        ZYN(JK) = ZYN(JWRK)
-      ENDDO
-    ENDDO
-    !    
-    !*      1.1.3   calcul des coordonnÊes gÊographiques des noeuds de la grille ISBA
-    !               -----------------------------------------------------------------
-    ALLOCATE(ZLAT(IMESHL))
-    ALLOCATE(ZLON(IMESHL))
-    CALL LATLON_CONF_PROJ(ZLAT0,ZLON0,ZRPK,ZBETA,ZLATOR,ZLONOR,ZXN,ZYN,ZLAT,ZLON)
-    DEALLOCATE(ZXN)
-    DEALLOCATE(ZYN)
-    !
-    !*      1.2 Gestion des coordonnees geographiques
-    !           -------------------------------------
-    !
-  ELSE IF(CGRID.EQ.'LONLAT REG') THEN
-    !
-    WRITE(ILUOUT,*) 'GRILLE LONLAT REG (application AMMA)' 
-    !
-    ALLOCATE(ZXI(NDIM_FULL))
-    ALLOCATE(ZYI(NDIM_FULL))
-    ZXI(:)=0.0
-    ZYI(:)=0.0
-    CALL GET_GRIDTYPE_LONLAT_REG(XGRID_PAR,PLONMIN=ZLONMIN,PLONMAX=ZLONMAX,             &
-                                 PLATMIN=ZLATMIN,PLATMAX=ZLATMAX,KLON=NIMAX,KLAT=NJMAX, &
-                                 KL=IL,PLON=ZXI,PLAT=ZYI)
-    !
-    IMESHL=(NIMAX+1)*(NJMAX+1)
-    !
-    ALLOCATE(ZLON(IMESHL))
-    ALLOCATE(ZLAT(IMESHL))
-    ALLOCATE(ZDXI(NDIM_FULL))
-    ALLOCATE(ZDYI(NDIM_FULL))
-    !
-    ZDXI(:)=(ZLONMAX-ZLONMIN)/(NIMAX-1)
-    ZDYI(:)=(ZLATMAX-ZLATMIN)/(NJMAX-1)
-    !
-    DO JJ=1,NIMAX
-      ZLON(JJ) = ZXI(JJ) - ZDXI(JJ)/2.
-    ENDDO
-    ZLON(NIMAX+1) = ZXI(NIMAX) + ZDXI(NIMAX)/2.
-    !
-    DO JJ=1,NJMAX
-      JWRK=(JJ-1)*(NIMAX+1)+1   ! indice sur la grille des noeuds
-      JI=(JJ-1)*NIMAX+1        ! indice sur la grille des mailles
-      ZLAT(JWRK) = ZYI(JI) - ZDYI(JI)/2.
-    ENDDO
-    !
-    JJ=((NJMAX+1)-1)*(NIMAX+1)+1  ! Indice sur la grille des noeuds
-    JI=(NJMAX-1)*NIMAX+1
-    ZLAT(JJ) = ZYI(JI) + ZDYI(JI)/2.
-    !
-    DEALLOCATE(ZDXI)
-    DEALLOCATE(ZDYI)
-    DEALLOCATE(ZXI)
-    DEALLOCATE(ZYI)
-    !
-    DO JJ=1,NIMAX+1
-      DO JI=2,NJMAX+1
-        JK = (JI-1)*(NIMAX+1)+JJ
-        ZLON(JK) = ZLON(JJ)
-      ENDDO
-    ENDDO
-    !
-    DO JI=1,NJMAX+1
-      DO JJ=2,NIMAX+1
-        JK=(JI-1)*(NIMAX+1)+JJ
-        JWRK=(JI-1)*(NIMAX+1)+1
-        ZLAT(JK)=ZLAT(JWRK)
-      ENDDO
-    ENDDO
-    !
-  ELSE
-    !       
-    WRITE(ILUOUT,*) 'ERREUR: TYPE DE GRILLE NON GERE PAR LE CODE'
-    CALL ABOR1_SFX("PGD_TOPD: TYPE DE GRILLE NON GERE PAR LE CODE")
-    !
-  ENDIF
-  !
-  !*      2.0  calcul des coordonnÊes lambert II Êtendu des noeuds de la grille ISBA
-  !            ----------------------------------------------------------------------
-  !
-  ALLOCATE(XXI(IMESHL))
-  ALLOCATE(XYI(IMESHL))
-  CALL XY_IGN(5,XXI,XYI,ZLAT,ZLON)
-  DEALLOCATE(ZLAT)
-  DEALLOCATE(ZLON)
-  !
-  !*      2.0     mask
-  !               ----
-  !***
-  CALL MAKE_MASK_TOPD_TO_ISBA(NDIM_FULL)
-  !
-  ALLOCATE(NNPIX(NDIM_FULL))
-  NNPIX(:) = NUNDEF
-  DO JJ=1,NDIM_FULL
-    NNPIX(JJ) = COUNT(NMASKT(:,:)==JJ)
-  ENDDO
-  !
-  CALL MAKE_MASK_ISBA_TO_TOPD(NDIM_FULL)
-  !
-  CALL WRITE_FILE_MASKTOPD(NDIM_FULL)
-  !
-  !*        3.0 Compute Mean slope over each ISBA_MESH
-  !            ----------------------------------------------------------------------
-  CALL TOPD_TO_ISBA_SLOPE(NDIM_FULL)
-  !
-  !*        4.0  Compute F and DC for each ISBA mesh
-  !            ----------------------------------------------------------------------
-  !
-  ALLOCATE(NNBV_IN_MESH(NDIM_FULL,NNCAT))
-  ALLOCATE(XBV_IN_MESH(NDIM_FULL,NNCAT))
-  ALLOCATE(XTOTBV_IN_MESH(NDIM_FULL))
-  !
-  XTOTBV_IN_MESH(:) = 0.0
-  !
-  DO JMESH=1,NDIM_FULL
-    XBV_IN_MESH(JMESH,:)=0.0
-    DO JCAT=1,NNCAT
-      NNBV_IN_MESH(JMESH,JCAT) = COUNT(NMASKI(JMESH,JCAT,:)/=NUNDEF)
-      XBV_IN_MESH(JMESH,JCAT) = REAL(NNBV_IN_MESH(JMESH,JCAT))*XDXT(JCAT)**2
-      XTOTBV_IN_MESH(JMESH) = XTOTBV_IN_MESH(JMESH) + XBV_IN_MESH(JMESH,JCAT)
-    ENDDO
-  ENDDO
-  !
-  ALLOCATE (ZF_PARAM(NDIM_FULL))
-  ALLOCATE (ZC_DEPTH_RATIO(NDIM_FULL))
-  !
-  ZF_PARAM(:) = 0.
-  ZC_DEPTH_RATIO(:) = 0.
-  DO JCAT=1,NNCAT
-    DO JMESH=1,NDIM_FULL
-      IF ( XTOTBV_IN_MESH(JMESH)/=0. ) THEN
-        ZF_PARAM(JMESH) = ZF_PARAM(JMESH) + XF_PARAM_BV(JCAT)*XBV_IN_MESH(JMESH,JCAT)/XTOTBV_IN_MESH(JMESH)
-        ZC_DEPTH_RATIO(JMESH) = ZC_DEPTH_RATIO(JMESH) + XC_DEPTH_RATIO_BV(JCAT)*XBV_IN_MESH(JMESH,JCAT)/XTOTBV_IN_MESH(JMESH)
-      ENDIF
-    ENDDO
-  ENDDO
-  !
-  WHERE (ZF_PARAM==0.)
-    ZF_PARAM=2.5
-    ZC_DEPTH_RATIO=1.
-  ENDWHERE
-  !
-  !write(*,*) 'f min max isba',MINVAL(ZF_PARAM),MAXVAL(ZF_PARAM)
-  !write(*,*) 'dc min max isba',MINVAL(ZC_DEPTH_RATIO),MAXVAL(ZC_DEPTH_RATIO)
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,'carte_f_dc.txt','FORMATTED',HACTION='WRITE')
-  DO JMESH=1,NDIM_FULL
-    WRITE(IUNIT,*) ZF_PARAM(JMESH),ZC_DEPTH_RATIO(JMESH)
-  ENDDO
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  !
-  DEALLOCATE(ZF_PARAM)
-  DEALLOCATE(ZC_DEPTH_RATIO)
-  !
-  WRITE(ILUOUT,*) 'Couplage avec TOPMODEL active'
-  !
-ELSE
-  !
-  WRITE(ILUOUT,*) 'Pas de couplage avec TOPMODEL'
-  !
-ENDIF
-!
-IF (LHOOK) CALL DR_HOOK('PGD_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE PGD_TOPD
diff --git a/src/SURFEX/prep_restart_coupl_topd.F90 b/src/SURFEX/prep_restart_coupl_topd.F90
deleted file mode 100644
index caf06a3bb23c47d7f02bf3445c77dead9e91b3e3..0000000000000000000000000000000000000000
--- a/src/SURFEX/prep_restart_coupl_topd.F90
+++ /dev/null
@@ -1,98 +0,0 @@
-!######
-SUBROUTINE PREP_RESTART_COUPL_TOPD(HPROGRAM,KI)
-!###################################################################
-!
-!!****  * PREP_RESTART_COUPL_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!!   
-!!    Write all files needed in case of restart of a simulation coupling SURFEX
-!!     and TOPODYN
-!!      
-!!    REFERENCE
-!!    ---------
-!!      
-!!    AUTHOR
-!!    ------
-!!	B. Vincendon    
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    07/06/11 
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, XQTOT, NNB_TOPD_STEP,&
-                                 XQB_RUN, XQB_DR
-USE MODD_COUPLING_TOPD, ONLY : XAS_NATURE,&
-                                 NNB_STP_RESTART, XWTOPT,&
-                                 XRUN_TOROUT, XDR_TOROUT
-!
-USE MODD_SURF_ATM_n,      ONLY:  NR_NATURE
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE MODI_WRITE_FILE_MAP
-USE MODI_UNPACK_SAME_RANK
-USE MODI_WRITE_FILE_ISBAMAP
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=6), INTENT(IN)         :: HPROGRAM ! program calling surf. schemes
-INTEGER,          INTENT(IN)         :: KI       ! Surfex grid dimension
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                        :: ILUOUT      ! unit of output listing file
-INTEGER                        :: IUNIT       ! unit of restart files
-INTEGER                        :: JSTP, JJ    ! loop control indexes
-REAL, DIMENSION(:),ALLOCATABLE :: ZAS         ! Saturated area fraction for each Isba meshes
- CHARACTER(LEN=30)              :: YVAR        ! name of results file
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('PREP_RESTART_COUPL_TOPD',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-! * 1. Write stock files
-!          
-WRITE(ILUOUT,*) 'Write STOCK file'
-!
- CALL OPEN_FILE('ASCII ',IUNIT,HFILE='stocks_sav.txt',HFORM='FORMATTED',HACTION='WRITE')
-DO JSTP = 1,NNB_STP_RESTART
-  WRITE(IUNIT,*)  XRUN_TOROUT(1:NNCAT,JSTP+NNB_TOPD_STEP), XDR_TOROUT(1:NNCAT,JSTP+NNB_TOPD_STEP)
-ENDDO
- CALL CLOSE_FILE('ASCII ',IUNIT)
-!  
-! * 2. Write pixels water content
-!
-WRITE(ILUOUT,*) 'Write pixels water content files'
-!
-YVAR = '_xwtop_sav.map'
- CALL WRITE_FILE_MAP(XWTOPT,YVAR)
-! 
-! * 3. Write Asat files
-! 
-WRITE(ILUOUT,*) 'Write Asat files'
-!
-ALLOCATE(ZAS(KI))
- CALL UNPACK_SAME_RANK(NR_NATURE,XAS_NATURE,ZAS)
-!
- CALL OPEN_FILE('ASCII ',IUNIT,HFILE='surfcont_sav.map',HFORM='FORMATTED',HACTION='WRITE')
- CALL WRITE_FILE_ISBAMAP(IUNIT,ZAS,KI)
- CALL CLOSE_FILE('ASCII ',IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('PREP_RESTART_COUPL_TOPD',1,ZHOOK_HANDLE)
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE PREP_RESTART_COUPL_TOPD
diff --git a/src/SURFEX/read_connex_file.F90 b/src/SURFEX/read_connex_file.F90
deleted file mode 100644
index 4561ed6a3ed4701d669b8f479890955e073afd9e..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_connex_file.F90
+++ /dev/null
@@ -1,108 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################     
-      SUBROUTINE READ_CONNEX_FILE(HPROGRAM,HFILE,HFORM,KNMC,PCONN,KLINE)
-!     #######################
-!
-!!****  *READ_CONNEX_FILE*  
-!!
-!!    PURPOSE
-!!    -------
-!     This routine aims at reading connexion file
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon    * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2006
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN, ONLY : NPMAX
-USE MODD_SURF_PAR,  ONLY : XUNDEF
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*),  INTENT(IN)  :: HPROGRAM    !
- CHARACTER(LEN=*),  INTENT(IN)  :: HFILE       ! File to be read
- CHARACTER(LEN=*),  INTENT(IN)  :: HFORM       ! Format of the file to be read
-INTEGER,           INTENT(IN)  :: KNMC       ! Number of pixels in the catchment
-REAL, DIMENSION(:,:),INTENT(OUT)   :: PCONN    ! pixels topographic slope/length flow
-INTEGER, DIMENSION(:),INTENT(OUT)  :: KLINE    ! second index of the pixel in the array PCONN
-!
-!*      0.2    declarations of local variables
-!
-!
-INTEGER                   :: JJ          ! loop control 
-INTEGER                   :: IUNIT       ! Unit of the files
-INTEGER                   :: ILUOUT      ! Unit of the files
-INTEGER                   :: IINDEX      ! index of the pixel in the topo domain
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_CONNEX_FILE',0,ZHOOK_HANDLE)
-!
-!*       0.2    preparing file openning
-!               ----------------------
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
- CALL OPEN_FILE(HPROGRAM,IUNIT,HFILE,HFORM,HACTION='READ')
-!
-WRITE(ILUOUT,*) 'Open ',HFILE,'debut'
-!
-DO JJ=1,7
-  READ(IUNIT,*) 
-ENDDO
-!
-DO JJ=1,KNMC
-  !
-  READ(IUNIT,'(f9.0,1x,f8.3,1x,2(f1.0,1x),8(f8.0,1x,f8.6,1x))',END=120) PCONN(JJ,:)
-  IINDEX = INT(PCONN(JJ,1))
-  KLINE(IINDEX) = JJ
-  !
-ENDDO
-!   
-120   CALL CLOSE_FILE(HPROGRAM,IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('READ_CONNEX_FILE',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_CONNEX_FILE
-
-
-
-
-
-
-
diff --git a/src/SURFEX/read_file_isbamap.F90 b/src/SURFEX/read_file_isbamap.F90
deleted file mode 100644
index 712a9d006820a7b4df97be5c77b71a92a90773b7..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_file_isbamap.F90
+++ /dev/null
@@ -1,101 +0,0 @@
-!-----------------------------------------------------------------
-!     ##########################
-      SUBROUTINE READ_FILE_ISBAMAP(KUNIT,PVAR,INI)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   25/01/2005
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-!
-USE MODD_TOPODYN
-USE MODD_SURF_ATM_GRID_n, ONLY : XGRID_PAR
-USE MODD_SURF_PAR, ONLY : XUNDEF
-!
-USE MODE_GRIDTYPE_CONF_PROJ
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN)             :: KUNIT  ! file unit
-REAL, DIMENSION(:), INTENT(OUT) :: PVAR   ! variable to write in the file
-INTEGER, INTENT(IN)             :: INI    ! Grid dimensions
-!
-!
-!*      0.2    declarations of local variables
-INTEGER                    :: JJ,JI
-INTEGER                    :: JINDEX ! reference number of the pixel
-REAL                       :: ZOUT
-REAL                       :: ZMAX,ZMIN
-REAL, DIMENSION(INI)       :: ZXI, ZYI   ! natural coordinates of ISBA grid (conformal projection)
-REAL, DIMENSION(INI)       :: ZDXI, ZDYI   ! Isba grid resolution in the conformal projection
-INTEGER                    :: IIMAX,IJMAX
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_FILE_ISBAMAP',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
- CALL GET_GRIDTYPE_CONF_PROJ(XGRID_PAR,PX=ZXI,PY=ZYI,KIMAX=IIMAX,KJMAX=IJMAX,PDX=ZDXI)
-!
-ZMAX=MAXVAL(PVAR)
-ZMIN=MINVAL(PVAR)
-ZOUT = XUNDEF
-!
-DO JJ=1,5
-  READ(KUNIT,*)
-ENDDO
-!
-READ(KUNIT,*) ZXI(1)
-READ(KUNIT,*) ZYI(1)
-READ(KUNIT,*) IIMAX
-READ(KUNIT,*) IJMAX
-READ(KUNIT,*) ZOUT
-READ(KUNIT,*) ZDXI(1)
-READ(KUNIT,*) ZMIN
-READ(KUNIT,*) ZMAX
-!
-DO JJ=1,IJMAX
-  DO JI=1,IIMAX
-    JINDEX=(JJ - 1) * IIMAX + JI
-    READ(KUNIT,*) PVAR(JINDEX)
-  ENDDO
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('READ_FILE_ISBAMAP',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_FILE_ISBAMAP
diff --git a/src/SURFEX/read_file_masktopd.F90 b/src/SURFEX/read_file_masktopd.F90
deleted file mode 100644
index 8425c8ecb3c129e7d597c4375b50013f5edd5abe..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_file_masktopd.F90
+++ /dev/null
@@ -1,121 +0,0 @@
-!-----------------------------------------------------------------
-!     ##########################
-      SUBROUTINE READ_FILE_MASKTOPD(KI)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2011
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-!
-USE MODD_TOPODYN
-USE MODD_COUPLING_TOPD, ONLY : NMASKI, NNPIX, NMASKT
-USE MODD_SURF_PAR,        ONLY : XUNDEF, NUNDEF
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-USE MODI_READ_TOPD_FILE
-!
-USE MODE_GRIDTYPE_CONF_PROJ
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN)             :: KI    ! Grid dimensions
-!
-!
-!*      0.2    declarations of local variables
-INTEGER                    :: JCAT,JMESH,JPIX
-INTEGER                    :: INUMPIX
-INTEGER :: IIMAX,IJMAX,IUNIT
- CHARACTER(LEN=50) :: YNAME
-REAL, DIMENSION(:),ALLOCATABLE    :: ZTOPD_READ !Topgraphic variable read
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_FILE_MASKTOPD',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
-ALLOCATE(ZTOPD_READ(NPMAX))
-!
-DO JCAT=1,NNCAT
-  !
-  YNAME=TRIM(CCAT(JCAT))//TRIM('.mask_topd')
-  !
-  CALL READ_TOPD_FILE('OFFLIN',YNAME,'FORMATTED',NNPT(JCAT),ZTOPD_READ)
-  !
-  DO JPIX=1,NPMAX
-    !
-    IF ( NLINE(JCAT,JPIX)/=0 ) THEN
-      IF  (ZTOPD_READ(JPIX)==XUNDEF ) THEN
-        NMASKT(JCAT,NLINE(JCAT,JPIX)) = NUNDEF
-      ELSE
-        NMASKT(JCAT,NLINE(JCAT,JPIX)) = FLOOR(ZTOPD_READ(JPIX))
-      ENDIF
-    ENDIF
-    !
-  ENDDO
-  !
-ENDDO
-!
-ALLOCATE(NNPIX(KI))
-NNPIX(:) = NUNDEF
-DO JMESH=1,KI
-  NNPIX(JMESH) = COUNT(NMASKT(:,:)==JMESH)
-ENDDO
-INUMPIX=MAXVAL(NNPIX)
-!
-ALLOCATE(NMASKI(KI,NNCAT,INUMPIX))
-NMASKI(:,:,:) = NUNDEF
-
-DO JCAT=1,NNCAT
-  !
-  YNAME=TRIM(CCAT(JCAT))//TRIM('.mask_surf')
-  CALL OPEN_FILE('ASCII ',IUNIT,YNAME,'FORMATTED','READ')
-  ! 
-  DO JMESH=1,KI
-    DO JPIX=1,NNPIX(JMESH)
-      READ(IUNIT,*) NMASKI(JMESH,JCAT,JPIX)
-    ENDDO
-  ENDDO
-  !
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  !
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('READ_FILE_MASKTOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_FILE_MASKTOPD
diff --git a/src/SURFEX/read_nam_pgd_topd.F90 b/src/SURFEX/read_nam_pgd_topd.F90
deleted file mode 100644
index efd433d63ef7a4392ac93030bb500f81617b4e0c..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_nam_pgd_topd.F90
+++ /dev/null
@@ -1,126 +0,0 @@
-!----------------------------------------------------------------------------!
-!     ##############################################################
-      SUBROUTINE READ_NAM_PGD_TOPD(HPROGRAM,OCOUPL_TOPD,HCAT,PF_PARAM_BV,PC_DEPTH_RATIO_BV)
-!     ##############################################################
-!
-!!**** *READ_NAM_TOPD_n* reads namelist NAM_TOPD
-!!
-!!    PURPOSE
-!!    -------
-!!    NAM_TOPD is a namelist used only for Topmodel coupling
-!!    It permits to define the different catchments studied.
-!!    This routine aims at reading and initialising those names.
-!!
-!!    METHOD
-!!    ------
-!!   
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!
-!!    B. Vincendon        Meteo-France
-!!
-!!    MODIFICATION
-!!    ------------
-!!
-!!    Original   11/2006
-!!
-!----------------------------------------------------------------------------
-!
-!*    0.     DECLARATION
-!            -----------
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_NAMELIST
-USE MODI_CLOSE_NAMELIST
-!
-USE MODD_TOPD_PAR, ONLY : JPCAT
-USE MODD_TOPODYN, ONLY : NNCAT
-!
-USE MODE_POS_SURF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*    0.1    Declaration of arguments
-!            ------------------------
-!
- CHARACTER(LEN=6),                   INTENT(IN)   :: HPROGRAM     ! Type of program
-LOGICAL, INTENT(OUT) :: OCOUPL_TOPD
- CHARACTER(LEN=15), DIMENSION(JPCAT),INTENT(OUT)  :: HCAT         ! Names of catchments         
-REAL, DIMENSION(JPCAT),INTENT(OUT)               :: PF_PARAM_BV
-REAL, DIMENSION(JPCAT),INTENT(OUT)               :: PC_DEPTH_RATIO_BV 
-!
-!*    0.2    Declaration of local variables
-!            ------------------------------
-!
- CHARACTER(LEN=15), DIMENSION(JPCAT) :: CCAT
-LOGICAL                           :: LCOUPL_TOPD
-REAL, DIMENSION(JPCAT)            :: XF_PARAM_BV
-REAL, DIMENSION(JPCAT)            :: XC_DEPTH_RATIO_BV
-!
-INTEGER                           :: ILUOUT    ! output listing logical unit
-INTEGER                           :: ILUNAM    ! namelist file logical unit
-LOGICAL                           :: GFOUND    ! flag when namelist is present
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-!*    0.3    Declaration of namelists
-!  
-!
-NAMELIST/NAM_PGD_TOPD/CCAT, LCOUPL_TOPD, XF_PARAM_BV, XC_DEPTH_RATIO_BV   
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_NAM_PGD_TOPD',0,ZHOOK_HANDLE)
-!
-!*    1.      Initializations of defaults
-!             ---------------------------
-!
-LCOUPL_TOPD = .FALSE.
-CCAT(:) = '   '
-XF_PARAM_BV(:) = 2.5
-XC_DEPTH_RATIO_BV(:) = 1.
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-!-------------------------------------------------------------------------------
-!
-!*    2.      Reading of namelist
-!             -------------------
-!
- CALL OPEN_NAMELIST(HPROGRAM,ILUNAM)
-!CALL OPEN_NAMELIST(HPROGRAM,'SURF  ',ILUNAM)
-!
- CALL POSNAM(ILUNAM,'NAM_PGD_TOPD',GFOUND,ILUOUT)
-IF (GFOUND) READ(UNIT=ILUNAM,NML=NAM_PGD_TOPD)
-!
- CALL CLOSE_NAMELIST(HPROGRAM,ILUNAM)
-!
-!         2.   Initialises number of catchments and time step variables
-!              -------------------------------------------------------
-!
-NNCAT=COUNT(CCAT(:)/='   ')
-!
-!-------------------------------------------------------------------------------
-!
-!*    3.      Fills output arguments
-!             ----------------------
-!
-OCOUPL_TOPD = LCOUPL_TOPD
-HCAT(1:NNCAT) = CCAT(1:NNCAT)
-PF_PARAM_BV(1:NNCAT) = XF_PARAM_BV(1:NNCAT)
-PC_DEPTH_RATIO_BV(1:NNCAT) = XC_DEPTH_RATIO_BV(1:NNCAT)
-!
-IF (LHOOK) CALL DR_HOOK('READ_NAM_PGD_TOPD',1,ZHOOK_HANDLE)
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE READ_NAM_PGD_TOPD
diff --git a/src/SURFEX/read_nam_topd.F90 b/src/SURFEX/read_nam_topd.F90
deleted file mode 100644
index 18893944a1e46aee509d75f718cbc1fd655accaa..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_nam_topd.F90
+++ /dev/null
@@ -1,158 +0,0 @@
-!----------------------------------------------------------------------------!
-!     ##############################################################
-       SUBROUTINE READ_NAM_TOPD(HPROGRAM,&
-                                OBUDGET_TOPD,KNB_TOPD,&
-                                OSTOCK_TOPD,&
-                                KNB_STOCK,KNB_RESTART,&
-                                PSPEEDR,PSPEEDG,PQINIT,PRTOP_D2)
-!     ##############################################################
-!
-!!**** *READ_NAM TOPD* reads namelist NAM_TOPD
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!    NAM_TOPD is a namelist used to define whether Topmodel coupling
-!!    is performed or not and the time step ratio between hydrological 
-!!    model and ISBA.
-!!    This routine aims at reading and initialising those variables.
-!!
-!!    METHOD
-!!    ------
-!!   
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!
-!!    B. Vincendon        Meteo-France
-!!
-!!    MODIFICATION
-!!    ------------
-!!
-!!    Original    11/2006
-!!
-!----------------------------------------------------------------------------
-!
-!*    0.     DECLARATION
-!            -----------
-!
-USE MODD_TOPD_PAR, ONLY : JPCAT
-USE MODD_TOPODYN, ONLY : NNCAT
-!
-USE MODE_POS_SURF
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_NAMELIST
-USE MODI_CLOSE_NAMELIST
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*    0.1    Declaration of arguments
-!            ------------------------
-!
- CHARACTER(LEN=6),                   INTENT(IN)   :: HPROGRAM     ! Type of program
-LOGICAL,                            INTENT(OUT)  :: OBUDGET_TOPD ! budget computation
-INTEGER,                            INTENT(OUT)  :: KNB_TOPD     ! Ratio between Topmodel time step and ISBA time step
-LOGICAL,                            INTENT(OUT)  :: OSTOCK_TOPD  ! T if use of stock from previous simulation
-INTEGER,                            INTENT(OUT)  :: KNB_STOCK    ! number of time step to read in previous simulation
-INTEGER,                            INTENT(OUT)  :: KNB_RESTART  ! number of time step to write for next simulation
-REAL, DIMENSION(JPCAT),INTENT(OUT)               :: PSPEEDR ! River speed
-REAL, DIMENSION(JPCAT),INTENT(OUT)               :: PSPEEDG ! Ground speed
-REAL, DIMENSION(JPCAT),INTENT(OUT)               :: PQINIT  ! Initial discharge at catchments outlet
-REAL, DIMENSION(JPCAT),INTENT(OUT)               :: PRTOP_D2
-!
-!*    0.2    Declaration of local variables
-!            ------------------------------
-!
-LOGICAL                           :: LBUDGET_TOPD
-LOGICAL                           :: LSTOCK_TOPD
-INTEGER                           :: NNB_TOPD
-INTEGER                           :: NFREQ_MAPS_WG
-INTEGER                           :: NFREQ_MAPS_ASAT
-INTEGER                           :: NNB_STP_STOCK
-INTEGER                           :: NNB_STP_RESTART
-REAL, DIMENSION(JPCAT)            :: XSPEEDR
-REAL, DIMENSION(JPCAT)            :: XSPEEDG
-REAL, DIMENSION(JPCAT)            :: XQINIT
-REAL, DIMENSION(JPCAT)            :: XRTOP_D2
-!
-INTEGER                           :: ILUOUT    ! output listing logical unit
-INTEGER                           :: ILUNAM    ! namelist file logical unit
-LOGICAL                           :: GFOUND    ! flag when namelist is present
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-!*    0.3    Declaration of namelists
-!  
-NAMELIST/NAM_TOPD/LBUDGET_TOPD, LSTOCK_TOPD, NNB_TOPD, &
-                  NFREQ_MAPS_WG, NFREQ_MAPS_ASAT,&
-                  NNB_STP_STOCK, NNB_STP_RESTART, &
-                  XSPEEDR, XSPEEDG, XQINIT, XRTOP_D2
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_NAM_TOPD',0,ZHOOK_HANDLE)
-!
-!*    1.      Initializations of defaults
-!             ---------------------------
-!
-LBUDGET_TOPD = .FALSE.
-LSTOCK_TOPD = .FALSE.
-NNB_TOPD = 1
-NFREQ_MAPS_WG = 0
-NFREQ_MAPS_ASAT = 0
-NNB_STP_STOCK = 1
-NNB_STP_RESTART = 1
-XSPEEDR(:) = 3.0 ! default value of river speed (adapted for Cevennes zone)
-XSPEEDG(:) = 0.3 ! default value of hillspeed
-XQINIT(:) = 0.
-XRTOP_D2(:) = 1.
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-!-------------------------------------------------------------------------------
-!
-!*    2.      Reading of namelist
-!             -------------------
-!
- CALL OPEN_NAMELIST(HPROGRAM,ILUNAM)
-!
- CALL POSNAM(ILUNAM,'NAM_TOPD',GFOUND,ILUOUT)
-IF (GFOUND) READ(UNIT=ILUNAM,NML=NAM_TOPD)
-!
- CALL CLOSE_NAMELIST(HPROGRAM,ILUNAM)
-!
-!-------------------------------------------------------------------------------
-!
-!*    3.      Fills output arguments
-!             ----------------------
-!
-OBUDGET_TOPD = LBUDGET_TOPD
-OSTOCK_TOPD = LSTOCK_TOPD
-KNB_TOPD = NNB_TOPD
-KNB_STOCK = NNB_STP_STOCK
-KNB_RESTART = NNB_STP_RESTART
-PSPEEDR(1:NNCAT) = XSPEEDR(1:NNCAT)
-PSPEEDG(1:NNCAT) = XSPEEDG(1:NNCAT)
-PQINIT(1:NNCAT) = XQINIT(1:NNCAT)
-PRTOP_D2(1:NNCAT) = XRTOP_D2(1:NNCAT)
-!
-WRITE(ILUOUT,*) 'NAM_TOPD:'
-WRITE(ILUOUT,*) 'LBUDGET ',LBUDGET_TOPD
-WRITE(ILUOUT,*) 'NNB_TOP',NNB_TOPD
-WRITE(ILUOUT,*) 'LSTOCK',LSTOCK_TOPD
-WRITE(ILUOUT,*) 'NNB_RESTART,NNB_STOCK',NNB_STP_RESTART,NNB_STP_STOCK
-!
-IF (LHOOK) CALL DR_HOOK('READ_NAM_TOPD',1,ZHOOK_HANDLE)
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE READ_NAM_TOPD
diff --git a/src/SURFEX/read_namelists_topd.F90 b/src/SURFEX/read_namelists_topd.F90
deleted file mode 100644
index 518b263e4067c61d9127e4fef1f0561aba9cd9f5..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_namelists_topd.F90
+++ /dev/null
@@ -1,39 +0,0 @@
-!     #########
-SUBROUTINE READ_NAMELISTS_TOPD(HPROGRAM)
-!     #######################################################
-!
-!--------------------------------------------------------------------------
-!
-USE MODI_READ_NAM_PGD_TOPD
-USE MODI_READ_NAM_TOPD
-!
-USE MODD_TOPODYN, ONLY : CCAT,XTOPD_STEP,NNB_TOPD_STEP,& 
-                           XSPEEDR,XQINIT,XRTOP_D2,XSPEEDG
-USE MODD_COUPLING_TOPD, ONLY : LCOUPL_TOPD, NNB_TOPD,&
-                                 LBUDGET_TOPD,LSTOCK_TOPD,&
-                                 NNB_STP_STOCK,NNB_STP_RESTART
-
-USE MODD_DUMMY_EXP_PROFILE,ONLY:XF_PARAM_BV,XC_DEPTH_RATIO_BV
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
- CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-!----------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('READ_NAMELISTS_TOPD',0,ZHOOK_HANDLE)
-!
- CALL READ_NAM_PGD_TOPD(HPROGRAM,LCOUPL_TOPD,CCAT,&
-                        XF_PARAM_BV,XC_DEPTH_RATIO_BV)
-IF (LCOUPL_TOPD) &
-  CALL READ_NAM_TOPD(HPROGRAM,LBUDGET_TOPD,NNB_TOPD,&
-                     LSTOCK_TOPD,NNB_STP_STOCK,NNB_STP_RESTART,&
-                     XSPEEDR,XSPEEDG,XQINIT,XRTOP_D2)
-!
-IF (LHOOK) CALL DR_HOOK('READ_NAMELISTS_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_NAMELISTS_TOPD
diff --git a/src/SURFEX/read_slope_file.F90 b/src/SURFEX/read_slope_file.F90
deleted file mode 100644
index 1177ad5998c9e908ef0b846c20a4f6043d71031f..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_slope_file.F90
+++ /dev/null
@@ -1,108 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################     
-      SUBROUTINE READ_SLOPE_FILE(HPROGRAM,HFILE,HFORM,KNMC,PTANB,PSLOP,PDAREA,PLAMBDA)
-!     #######################
-!
-!!****  *READ_SLOPE_FILE*  
-!!
-!!    PURPOSE
-!!    -------
-!     This routine aims at reading topographic files
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon    * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2006
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE MODD_TOPODYN, ONLY : NPMAX
-USE MODD_SURF_PAR,  ONLY : XUNDEF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*),  INTENT(IN)  :: HPROGRAM    !
- CHARACTER(LEN=*),  INTENT(IN)  :: HFILE       ! File to be read
- CHARACTER(LEN=*),  INTENT(IN)  :: HFORM       ! Format of the file to be read
-INTEGER,           INTENT(IN)  :: KNMC        ! Number of pixels in the catchment
-REAL, DIMENSION(:),  INTENT(OUT)   :: PTANB    ! pixels topographic slope(tan(beta)
-REAL, DIMENSION(:),  INTENT(OUT)   :: PSLOP   ! pixels topographic slope/length flow
-REAL, DIMENSION(:),  INTENT(OUT)   :: PDAREA  ! drainage area (aire drainee)
-REAL, DIMENSION(:),  INTENT(OUT)   :: PLAMBDA ! pure topographic index
-!
-!*      0.2    declarations of local variables
-!
-!
-INTEGER                   :: JJ ! loop control 
-INTEGER                   :: IWRK        ! work variable
-INTEGER                   :: IUNIT       ! Unit of the files
-INTEGER                   :: ILUOUT      ! Unit of the files
-!
-REAL                      :: ZWRK        ! work variable
-REAL, DIMENSION(KNMC)     :: ZDAREA      ! drainage area (aire drainee)
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_SLOPE_FILE',0,ZHOOK_HANDLE)
-!
-!*       0.2    preparing file openning
-!               ----------------------
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
- CALL OPEN_FILE(HPROGRAM,IUNIT,HFILE,HFORM,HACTION='READ')
-!
-READ(IUNIT,*) 
-!
-DO JJ=1,KNMC
-  !
-   READ(IUNIT,*,END=110) IWRK, PTANB(JJ), PSLOP(JJ), ZWRK, PDAREA(JJ)
-  PLAMBDA(JJ) = LOG(PDAREA(JJ)/PSLOP(JJ))
-  !
-ENDDO
-!
-110 CALL CLOSE_FILE(HPROGRAM,IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('READ_SLOPE_FILE',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_SLOPE_FILE
-
-
-
-
-
-
-
diff --git a/src/SURFEX/read_topd_file.F90 b/src/SURFEX/read_topd_file.F90
deleted file mode 100644
index e1f84ec009410d87fd09b8a6e15b98bb3452d658..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_topd_file.F90
+++ /dev/null
@@ -1,105 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE READ_TOPD_FILE(HPROGRAM,HFILE,HFORM,KNPT,PTOPD_READ)
-!     #######################
-!
-!!****  *READ_TOPD_FILE*  
-!!
-!!    PURPOSE
-!!    -------
-!     This routine aims at reading topographic files
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon    * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2006
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE MODD_TOPODYN, ONLY : NPMAX
-USE MODD_SURF_PAR,  ONLY : XUNDEF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*),  INTENT(IN)  :: HPROGRAM    !
- CHARACTER(LEN=*),  INTENT(IN)  :: HFILE       ! File to be read
- CHARACTER(LEN=*),  INTENT(IN)  :: HFORM       ! Format of the file to be read
-INTEGER,           INTENT(IN)  :: KNPT        ! Number of points in the catchment
-REAL, DIMENSION(:), INTENT(OUT) :: PTOPD_READ ! Topographic parameter read on file
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                   :: JJ,JI  ! loop control
-INTEGER                   :: IUNIT  ! Unit of the files
-INTEGER                   :: ILUOUT ! Unit of the files
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_TOPD_FILE',0,ZHOOK_HANDLE)
-!
-!*       0.2    preparing file openning
-!               ----------------------
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
- CALL OPEN_FILE(HPROGRAM,IUNIT,HFILE,HFORM,HACTION='READ')
-!PTOPD_READ(:)=XUNDEF
-!
-DO JJ=1,13
-  READ(IUNIT,*) 
-ENDDO
-!
-DO JI=1,KNPT
-  !
-  READ(IUNIT,*,END=110) PTOPD_READ(JI)
-  !
-  IF (PTOPD_READ(JI).LT.0.0) PTOPD_READ(JI) = XUNDEF
-  !
-ENDDO
-!
-110 CALL CLOSE_FILE(HPROGRAM,IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('READ_TOPD_FILE',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_TOPD_FILE
-
-
-
-
-
-
-
diff --git a/src/SURFEX/read_topd_header_connex.F90 b/src/SURFEX/read_topd_header_connex.F90
deleted file mode 100644
index ea19f2eed236b8ae91c8ea8d7301a3f205b12917..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_topd_header_connex.F90
+++ /dev/null
@@ -1,100 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE READ_TOPD_HEADER_CONNEX(HPROGRAM,HFILE,HFORM,KNMC)
-!     #######################
-!
-!!****  *READ_TOPD_HEADER*  
-!!
-!!    PURPOSE
-!!    -------
-!     This routine aims at reading topographic files
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon    * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2006
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE MODD_TOPODYN, ONLY : NPMAX
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*),  INTENT(IN)  :: HPROGRAM    !
- CHARACTER(LEN=*),  INTENT(IN)  :: HFILE       ! File to be read
- CHARACTER(LEN=*),  INTENT(IN)  :: HFORM       ! Format of the file to be read
-INTEGER,           INTENT(OUT) :: KNMC     ! number of pixels in a catchment 
-!
-!*      0.2    declarations of local variables
-!
-!
-INTEGER                   :: JJ ! loop control 
-INTEGER                   :: IUNIT       ! Unit of the files
-INTEGER                   :: ILUOUT      ! Unit of the files
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_TOPD_HEADER_CONNEX',0,ZHOOK_HANDLE)
-!
-!*       0.2    preparing file openning
-!               ----------------------
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-WRITE(ILUOUT,*) 'Open ',HFILE,'header'
-!
- CALL OPEN_FILE(HPROGRAM,IUNIT,HFILE,HFORM,HACTION='READ')
-!
-READ(IUNIT,*)
-READ(IUNIT,*) KNMC
-!
-DO JJ=1,5
-  READ(IUNIT,*) 
-ENDDO
-!
- CALL CLOSE_FILE(HPROGRAM,IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('READ_TOPD_HEADER_CONNEX',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_TOPD_HEADER_CONNEX
-
-
-
-
-
-
- 
diff --git a/src/SURFEX/read_topd_header_dtm.F90 b/src/SURFEX/read_topd_header_dtm.F90
deleted file mode 100644
index 5e192db15c07df2630eecacc9980b03b81b3765e..0000000000000000000000000000000000000000
--- a/src/SURFEX/read_topd_header_dtm.F90
+++ /dev/null
@@ -1,111 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE READ_TOPD_HEADER_DTM(HPROGRAM,HFILE,HFORM,PX0,PY0,KNXC,KNYC,PNUL,PDXT)
-!     #######################
-!
-!!****  *READ_TOPD_HEADER*  
-!!
-!!    PURPOSE
-!!    -------
-!     This routine aims at reading topographic files header 
-!     for a given file (then for a cathment)
-!
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon    * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2006
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE MODD_TOPODYN, ONLY : NPMAX
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*),  INTENT(IN)  :: HPROGRAM   !
- CHARACTER(LEN=*),  INTENT(IN)  :: HFILE      ! File to be read
- CHARACTER(LEN=*),  INTENT(IN)  :: HFORM      ! Format of the file to be read
-REAL,              INTENT(OUT) :: PX0        ! abcissa of bottom-left pixel   
-REAL,              INTENT(OUT) :: PY0        ! ordinate of bottom-left pixel   
-INTEGER,           INTENT(OUT) :: KNXC       ! number of topographixc grid points along abcissa axis   
-INTEGER,           INTENT(OUT) :: KNYC       ! number of topographixc grid points along ordinate axis   
-REAL,              INTENT(OUT) :: PNUL       ! undifined value in topographic files
-REAL,              INTENT(OUT) :: PDXT       ! catchment rid mesh size   
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                   :: JJ          ! loop control 
-INTEGER                   :: IUNIT       ! Unit of the files
-INTEGER                   :: ILUOUT      ! Unit of the files
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('READ_TOPD_HEADER_DTM',0,ZHOOK_HANDLE)
-!
-!*       0.2    preparing file openning
-!               ----------------------
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-WRITE(ILUOUT,*) 'Open ',HFILE,'header'
-!
- CALL OPEN_FILE(HPROGRAM,IUNIT,HFILE,HFORM,HACTION='READ')
-!
-DO JJ=1,5
-  READ(IUNIT,*) 
-ENDDO
-!
-READ(IUNIT,*) PX0
-READ(IUNIT,*) PY0
-READ(IUNIT,*) KNXC
-READ(IUNIT,*) KNYC
-READ(IUNIT,*) PNUL
-READ(IUNIT,*) PDXT
-READ(IUNIT,*)
-READ(IUNIT,*)
-!
- CALL CLOSE_FILE(HPROGRAM,IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('READ_TOPD_HEADER_DTM',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE READ_TOPD_HEADER_DTM
-
-
-
-
-
-
- 
diff --git a/src/SURFEX/recharge_surf_topd.F90 b/src/SURFEX/recharge_surf_topd.F90
deleted file mode 100644
index e28d2db7227c31895bc80fdf067c0c7bf7e21086..0000000000000000000000000000000000000000
--- a/src/SURFEX/recharge_surf_topd.F90
+++ /dev/null
@@ -1,230 +0,0 @@
-!-----------------------------------------------------------
-!     #######################
-      SUBROUTINE RECHARGE_SURF_TOPD(PHI,PHT,KI)
-!     #######################
-!
-!!****  *RECHARGE_SURF_TOPD*  
-!!
-!!    PURPOSE
-!!    -------
-!
-!     
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   12/2003
-!! 
-!!    WARNING
-!!    ----------------
-!!     on considčre que le seuil pour les deficits reste le niveau Wfc
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-USE MODD_COUPLING_TOPD, ONLY: NMASKI, XWFCTOPT, XDMAXFC, XWTOPT,&
-                                XDTOPT, XWSTOPT, NNPIX
-USE MODD_TOPODYN,       ONLY: NNCAT, XDMAXT
-!
-USE MODD_SURF_PAR,        ONLY: XUNDEF,NUNDEF
-!
-USE MODI_ABOR1_SFX
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-!
-INTEGER, INTENT(IN) :: KI    ! Grid dimensions
-REAL, DIMENSION(:), INTENT(INOUT)   :: PHI   ! water content variation since last time step from ISBA (m)
-REAL, DIMENSION(:,:), INTENT(OUT)   :: PHT   ! water content variation to provide to TOPODYN to be distributed (m) 
-!
-!*      0.2    declarations of local variables
-!
-!
-LOGICAL, DIMENSION(NNCAT,SIZE(NMASKI,3)) :: GTEST
-INTEGER            :: J1,J2,J3,J4     ! loop control 
-INTEGER            :: INBSAT, INBALL
-!
-REAL                           :: ZREST            ! m
-REAL                           :: ZWNEW            ! m3/m3
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('RECHARGE_SURF_TOPD',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
-!*       1.     ISBA => TOPODYN-LAT
-!               -------------------
-!
-PHT(:,:)=0.
-!
-DO J3 = 1,KI
-  !
-  !The water content is lower than the previous one : this case is dealed with first to fasten the computation time.
-  IF (PHI(J3) <= 0.0) THEN
-    !
-    DO J1 = 1,NNCAT
-      !
-      J4 = 1
-      J2 = NMASKI(J3,J1,J4)
-      !
-      DO WHILE (J2 /= NUNDEF .AND. J4<=SIZE(NMASKI,3) )
-        !
-        IF ( XWFCTOPT(J1,J2) /= XUNDEF ) THEN
-          !
-          ZWNEW = XWTOPT(J1,J2) + PHI(J3) / XDTOPT(J1,J2)
-          !
-          IF ( ZWNEW >= XWFCTOPT(J1,J2) ) THEN
-            !
-            ! on reste au-dessus de la capacite au champ, malgre l'assechement
-            XDMAXT(J1,J2) = XDMAXFC(J1,J2) 
-            PHT(J1,J2) = (ZWNEW - XWFCTOPT(J1,J2)) * XDTOPT(J1,J2)
-            !
-          ELSE ! on passe au-dessous de la capacite au champ
-            !
-            XDMAXT(J1,J2) = (XWSTOPT(J1,J2) - ZWNEW) * XDTOPT(J1,J2)
-            PHT(J1,J2) = 0.0
-            !
-          ENDIF
-          !
-          J4 = J4+1
-          IF ( J4<=SIZE(NMASKI,3) ) J2 = NMASKI(J3,J1,J4)
-          !
-        ELSE ! pixel non défini dans Isba
-          !
-          XDMAXT(J1,J2) = 0.0
-          PHT = 0.0
-          !
-        ENDIF
-        !
-      ENDDO
-      !
-    ENDDO
-    !
-  ELSE ! recharge > 0.0
-    !
-    ZREST=1.
-    GTEST(:,:)=.TRUE.
-    !
-    DO WHILE ( ZREST>0.0 )
-      !
-      ZREST=0.0
-      !
-      DO J1=1,NNCAT
-        !
-        J4=1
-        J2=NMASKI(J3,J1,J4)
-        !
-        DO WHILE ( J2/=NUNDEF .AND. J4<=SIZE(NMASKI,3) )
-          !
-          IF ( GTEST(J1,J4) .AND. XWFCTOPT(J1,J2)/=XUNDEF ) THEN
-            !
-            ZWNEW = XWTOPT(J1,J2) + PHI(J3) / XDTOPT(J1,J2)
-            !
-            IF ( XWTOPT(J1,J2) == XWSTOPT(J1,J2) ) THEN ! pixel déjā saturé
-              !
-              XDMAXT(J1,J2) = 0.0
-              PHT(J1,J2) = 0.0
-              ZREST = ZREST + PHI(J3)
-              GTEST(J1,J4) = .FALSE.
-              !
-            ELSE IF ( ( XWSTOPT(J1,J2) - XWTOPT(J1,J2) ) * XDTOPT(J1,J2) <= PHI(J3) ) THEN
-              !
-              ! pixel va se saturer
-              XDMAXT(J1,J2) = XDMAXFC(J1,J2)
-              PHT(J1,J2) = ( XWSTOPT(J1,J2) - XWFCTOPT(J1,J2) ) * XDTOPT(J1,J2)
-              ZREST = ZREST + PHI(J3) - PHT(J1,J2)
-              GTEST(J1,J4)=.FALSE.
-              !
-            ELSE IF ( XWTOPT(J1,J2) < XWFCTOPT(J1,J2) ) THEN 
-              !
-              ! en dessous de la capacité au champ avant d'ajouter la recharge
-              IF ( (XWTOPT(J1,J2) + PHI(J3)/XDTOPT(J1,J2)) <= XWFCTOPT(J1,J2) ) THEN 
-                !
-                ! en dessous de la capacité au champ avec la recharge 
-                XDMAXT(J1,J2) = ( XWSTOPT(J1,J2) - ZWNEW ) * XDTOPT(J1,J2)
-                PHT(J1,J2) = 0.0
-                !
-              ELSE ! au-dessus de la capacité au champ avec la recharge
-                !
-                XDMAXT(J1,J2) = XDMAXFC(J1,J2)
-                PHT(J1,J2) = ( ZWNEW - XWFCTOPT(J1,J2) ) * XDTOPT(J1,J2)
-                !
-              ENDIF
-              !
-            ELSE ! au-dessus de la capacité au champ avant d'ajouter la recharge
-              !
-              XDMAXT(J1,J2) = XDMAXFC(J1,J2)
-              PHT(J1,J2) = ( ZWNEW - XWFCTOPT(J1,J2) ) * XDTOPT(J1,J2)
-              !
-            ENDIF
-            !
-          ELSE IF ( XWFCTOPT(J1,J2)==XUNDEF ) THEN! pixel non défini dans Isba
-            !
-            XDMAXT(J1,J2) = 0.0
-            PHT = 0.0
-            !
-          ENDIF
-          !
-          J4 = J4+1
-          IF ( J4<=SIZE(NMASKI,3) ) J2 = NMASKI(J3,J1,J4)
-          !
-        ENDDO
-        !
-      ENDDO
-      !
-      IF ( ZREST/=0.0 ) THEN
-        !
-        INBSAT=COUNT(.NOT.GTEST) !nb de pixels saturés avec ou sans la recharge
-        !
-        IF ( INBSAT == NNPIX(J3) ) THEN
-          !
-          IF (NNPIX(J3) > 400 ) THEN
-            WRITE(*,*) 'MAILLE NUM=',J3, 'nb pix tot=',NNPIX(J3)
-            CALL ABOR1_SFX("RECHARGE_SURF_TOPD: TOO MANY PIXELS SATURATED ")
-          ELSE
-            ZREST=0.0
-          ENDIF
-          !
-        ELSE
-          !
-          PHI(J3) = PHI(J3) + ( ZREST / (NNPIX(J3) - INBSAT) ) ! nouvelle recharge ā distribuer
-          !
-        ENDIF
-      ENDIF
-      !
-    ENDDO
-    !
-  ENDIF
-  !
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('RECHARGE_SURF_TOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE RECHARGE_SURF_TOPD
diff --git a/src/SURFEX/restart_coupl_topd.F90 b/src/SURFEX/restart_coupl_topd.F90
deleted file mode 100644
index 0eccccc0654015a6686558d30adb0396a0ba2c35..0000000000000000000000000000000000000000
--- a/src/SURFEX/restart_coupl_topd.F90
+++ /dev/null
@@ -1,180 +0,0 @@
-!######
-SUBROUTINE RESTART_COUPL_TOPD(HPROGRAM,KI)
-!###################################################################
-!
-!!****  *RESTART_COUPL_TOPDn*  
-!!
-!!    PURPOSE
-!!    -------
-!!  Read all files needed in case of restart 
-!!      
-!!    REFERENCE
-!!    ---------
-!!      
-!!    AUTHOR
-!!    ------
-!!	B. Vincendon    
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    07/06/11 
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_SURF_PAR,        ONLY : XUNDEF,NUNDEF
-USE MODD_SURF_ATM_n,      ONLY:  NR_NATURE
-USE MODD_ISBA_PAR,        ONLY : XWGMIN
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, CCAT, NNPT, NLINE, NNMC, NPMAX,&
-                               NNB_TOPD_STEP
-
-USE MODD_COUPLING_TOPD, ONLY : XAS_NATURE, &
-                                 NNB_STP_STOCK,NNB_STP_RESTART,XWTOPT,&
-                                 XRUN_TOROUT,XDR_TOROUT
-!
-USE MODI_READ_TOPD_FILE
-USE MODI_READ_FILE_ISBAMAP
-!
-USE MODI_GET_LUOUT
-USE MODI_ABOR1_SFX
-USE MODI_UNPACK_SAME_RANK
-USE MODI_PACK_SAME_RANK
-!
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=6), INTENT(IN)         :: HPROGRAM ! program calling surf. schemes
-INTEGER,          INTENT(IN)         :: KI       ! Surfex grid dimension
-!
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                                     :: ILUOUT   ! unit of output listing file
-INTEGER                                     :: IUNIT    ! unit of restart files
-INTEGER                                     :: JSTP,JCAT,JPIX! loop control indexes
-REAL, DIMENSION(:),ALLOCATABLE              :: ZAS      ! Saturated area fraction for each Isba meshes
-REAL, DIMENSION(:),ALLOCATABLE              :: ZWTOPT   ! Initial water content in case of restart
- CHARACTER(LEN=50), DIMENSION(:),ALLOCATABLE :: YFILETOP ! File names
-LOGICAL                                     :: LSTOCK, LWG, LASAT
-REAL                                        :: ZCORR_STOCK ! used to avoid to lose stock
-REAL                                        :: ZCNT_UNDEF,ZSUM1,ZSUM2 ! used to correct budget
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('RESTART_COUPL_TOPD',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-! * 1. Read stock files
-!          
-WRITE(*,*) 'Read STOCK file ',NNB_STP_STOCK
-NNB_STP_STOCK = MIN(NNB_STP_STOCK, NNB_TOPD_STEP + NNB_STP_RESTART)
-!
-INQUIRE(FILE='stocks_init.txt', EXIST=LSTOCK)
-INQUIRE(FILE='surfcont_init.map', EXIST=LASAT)
-!
-IF (.NOT.LSTOCK) THEN
-  WRITE(ILUOUT,*) 'You asked to run in restart mode but stock file is missing'
-  CALL ABOR1_SFX("RESTART_COUPL_TOPD_n: stock file is missing")
-ELSEIF (.NOT.LASAT) THEN
-  WRITE(ILUOUT,*) 'You asked to run in restart mode but contributive area file is missing'
-  CALL ABOR1_SFX("RESTART_COUPL_TOPD_n: contributive area file is missing")
-ELSE
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,'stocks_init.txt','FORMATTED',HACTION='READ ')
-  DO JSTP=1,NNB_STP_STOCK
-    READ(IUNIT,*)  XRUN_TOROUT(1:NNCAT,JSTP),XDR_TOROUT(1:NNCAT,JSTP)
-  ENDDO
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  !  
-  ! * 2. Read pixels water content
-  !
-  DO JCAT=1,NNCAT
-    !
-    IF (XRUN_TOROUT(JCAT,NNB_STP_STOCK)/=0.) THEN
-      !
-      ZCORR_STOCK = XRUN_TOROUT(JCAT,NNB_STP_STOCK) - XRUN_TOROUT(JCAT,NNB_STP_STOCK-1)
-      DO JSTP = NNB_STP_STOCK+1,NNB_TOPD_STEP
-        XRUN_TOROUT(JCAT,JSTP) = MAX(0.,XRUN_TOROUT(JCAT,JSTP-1)+ZCORR_STOCK)
-      ENDDO
-      !
-    ENDIF
-    !
-    IF (XDR_TOROUT(JCAT,NNB_STP_STOCK)/=0.) THEN
-      !
-      ZCORR_STOCK = XDR_TOROUT(JCAT,NNB_STP_STOCK) - XDR_TOROUT(JCAT,NNB_STP_STOCK-1)
-      DO JSTP = NNB_STP_STOCK+1,NNB_TOPD_STEP
-        XDR_TOROUT(JCAT,JSTP) = MAX(0.,XDR_TOROUT(JCAT,JSTP-1)+ZCORR_STOCK)
-      ENDDO
-      !
-    ENDIF
-    !
-  ENDDO
-  !
-  WRITE(*,*) 'Write pixels water content files'
-  !
-  ALLOCATE(ZWTOPT(NPMAX))
-  ALLOCATE(YFILETOP(NNCAT))
-  !
-  DO JCAT=1,NNCAT
-    !
-    YFILETOP(JCAT)=TRIM(CCAT(JCAT))//'_xwtop_init.map'
-    INQUIRE(FILE=YFILETOP(JCAT), EXIST=LWG)
-    IF (.NOT.LWG) THEN
-      !
-      WRITE(ILUOUT,*) 'You asked to run in restart mode but pixels water content file is missing'
-      WRITE(ILUOUT,*) 'for catchment : ',CCAT(JCAT)
-      CALL ABOR1_SFX("RESTART_COUPL_TOPD_n: pixels water content file is missing")
-      !
-    ELSE
-      !
-      ZSUM1=SUM(XWTOPT(JCAT,:),MASK=XWTOPT(JCAT,:)/=XUNDEF)
-      !
-      CALL READ_TOPD_FILE('ASCII ',YFILETOP(JCAT),'FORMATTED',NNPT(JCAT),ZWTOPT)
-      !
-      DO JPIX=1,SIZE(NLINE(JCAT,:))
-        IF ( NLINE(JCAT,JPIX)/=0 .AND. NLINE(JCAT,JPIX)/=XUNDEF ) THEN
-          IF (ZWTOPT(JPIX) /= XUNDEF) XWTOPT(JCAT,NLINE(JCAT,JPIX)) = ZWTOPT(JPIX)
-        ENDIF
-      ENDDO
-      !
-      ZSUM2=SUM(XWTOPT(JCAT,:),MASK=XWTOPT(JCAT,:)<XUNDEF)
-      !
-      IF ( ABS(ZSUM2-ZSUM1)>100. ) THEN
-        !
-        ZCNT_UNDEF = COUNT(XWTOPT(JCAT,NLINE(JCAT,:))/=XUNDEF.AND. NLINE(JCAT,:)/=0)
-        IF (ZCNT_UNDEF/=0.) THEN
-          WHERE ( XWTOPT(JCAT,NLINE(JCAT,:))/=XUNDEF .AND. NLINE(JCAT,:)/=0 )
-            XWTOPT(JCAT,NLINE(JCAT,:)) = XWTOPT(JCAT,NLINE(JCAT,:)) - ((ZSUM2-ZSUM1)/ZCNT_UNDEF)
-          ENDWHERE
-        ENDIF
-        ZSUM2=SUM(XWTOPT(JCAT,:),MASK=XWTOPT(JCAT,:)<XUNDEF)
-        !
-      ENDIF
-      !
-    ENDIF
-    !
-  ENDDO
-  ! 
-  ! * 3. Read Asat files
-  ! 
-  WRITE(*,*) 'Write Asat files'
-  ALLOCATE(ZAS(KI))
-  CALL OPEN_FILE('ASCII ',IUNIT,'surfcont_init.map','FORMATTED',HACTION='READ ')
-  CALL READ_FILE_ISBAMAP(IUNIT,ZAS,KI)
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  CALL PACK_SAME_RANK(NR_NATURE,ZAS,XAS_NATURE)
-  !
-ENDIF
-!
-IF (LHOOK) CALL DR_HOOK('RESTART_COUPL_TOPD',1,ZHOOK_HANDLE)
-!-------------------------------------------------------------------------------
-END SUBROUTINE RESTART_COUPL_TOPD
diff --git a/src/SURFEX/rout_data_isba.F90 b/src/SURFEX/rout_data_isba.F90
deleted file mode 100644
index 2bb8abc3c97d9cff4f36e20c6a12a3a57d6e2adc..0000000000000000000000000000000000000000
--- a/src/SURFEX/rout_data_isba.F90
+++ /dev/null
@@ -1,146 +0,0 @@
-!-----------------------------------------------------------------
-!     #####################
-      SUBROUTINE ROUT_DATA_ISBA(HPROGRAM,KI,KSTEP)
-!     #####################
-!
-!!****  *ROUT_DATA_ISBA*  
-!!
-!!    PURPOSE
-!!    -------
-!
-!    Routes runoff and drainage of ISBA with coupling with Topmodel
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon	*  Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   15/06/2007
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODI_GET_LUOUT
-USE MODI_UNPACK_SAME_RANK
-USE MODI_DIAG_ISBA_TO_ROUT
-USE MODI_ISBA_TO_TOPD
-USE MODI_ROUTING
-!
-USE MODD_TOPODYN,        ONLY : NNCAT, NMESHT, NNMC
-USE MODD_COUPLING_TOPD,  ONLY : NMASKT, XRUNOFF_TOP, XATOP, NNPIX,&
-                                  XAVG_RUNOFFCM, XAVG_DRAINCM
-!
-USE MODD_DIAG_EVAP_ISBA_n, ONLY : XAVG_DRAINC
-USE MODD_SURF_ATM_n,       ONLY : NR_NATURE
-USE MODD_SURF_PAR,         ONLY : XUNDEF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
-INTEGER, INTENT(IN)          :: KI     ! Grid dimensions
-INTEGER, INTENT(IN)          :: KSTEP  ! current time step 
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                       :: JJ,JI  ! loop control 
-INTEGER                       :: IUNIT       ! unit number of results files
-INTEGER                       :: ILUOUT      ! unit number of listing file
- CHARACTER(LEN=30)             :: YVAR
-REAL, DIMENSION(KI)           :: ZRUNOFFC_FULL  ! Cumulated runoff from isba on the full domain (kg/m2)
-REAL, DIMENSION(KI)           :: ZRUNOFFC_FULLM ! Cumulated runoff from isba on the full domain (kg/m2) at t-dt
-REAL, DIMENSION(KI)           :: ZRUNOFF_ISBA   ! Runoff from Isba (kg/m2)
-REAL, DIMENSION(KI)           :: ZDRAINC_FULL   ! Cumulated drainage from Isba on the full domain (kg/m2)
-REAL, DIMENSION(KI)           :: ZDRAINC_FULLM  ! Cumulated drainage from Isba on the full domain (kg/m2) at t-dt
-REAL, DIMENSION(KI)           :: ZDRAIN_ISBA    ! Drainage from Isba (m3/s)
-REAL, DIMENSION(NNCAT,NMESHT) :: ZRUNOFF_TOPD   ! Runoff on the Topodyn grid (m3/s)
-REAL, DIMENSION(NNCAT,NMESHT) :: ZDRAIN_TOPD    ! Drainage from Isba on Topodyn grid (m3/s)
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('ROUT_DATA_ISBA',0,ZHOOK_HANDLE)
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-!
-ZRUNOFFC_FULL (:) = 0.
-ZRUNOFFC_FULLM(:) = 0.
-ZRUNOFF_ISBA  (:) = 0.
-ZRUNOFF_TOPD(:,:) = 0.
-ZDRAINC_FULL  (:) = 0.
-ZDRAINC_FULLM (:) = 0.
-ZDRAIN_ISBA   (:) = 0.
-ZDRAIN_TOPD (:,:) = 0.
-!
-!    Runoff on TOPODYN grid
-!   ---------------------------------------
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XRUNOFF_TOP,ZRUNOFFC_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_RUNOFFCM,ZRUNOFFC_FULLM)
-!
- CALL DIAG_ISBA_TO_ROUT(ZRUNOFFC_FULL,ZRUNOFFC_FULLM,ZRUNOFF_ISBA)
-!
-XAVG_RUNOFFCM(:) = XRUNOFF_TOP(:)
-ZRUNOFF_TOPD(:,:) = 0.0
-!
- CALL ISBA_TO_TOPD(ZRUNOFF_ISBA,ZRUNOFF_TOPD)
-!
-DO JJ=1,NNCAT
-  DO JI=1,NNMC(JJ)
-    ZRUNOFF_TOPD(JJ,JI) = ZRUNOFF_TOPD(JJ,JI) / NNPIX(NMASKT(JJ,JI))
-  ENDDO
-ENDDO
-!
-!    Drainage treatment
-!    ----------------------------------------
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINC*XATOP,ZDRAINC_FULL)
- CALL UNPACK_SAME_RANK(NR_NATURE,XAVG_DRAINCM*XATOP,ZDRAINC_FULLM)
-!
- CALL DIAG_ISBA_TO_ROUT(ZDRAINC_FULL,ZDRAINC_FULLM,ZDRAIN_ISBA)
-!
-XAVG_DRAINCM(:)  = XAVG_DRAINC(:)
-ZDRAIN_TOPD(:,:) = 0.0
-!
- CALL ISBA_TO_TOPD(ZDRAIN_ISBA,ZDRAIN_TOPD)
-!
-DO JJ=1,NNCAT
-  DO JI=1,NNMC(JJ)
-    ZDRAIN_TOPD(JJ,JI) = ZDRAIN_TOPD(JJ,JI) / NNPIX(NMASKT(JJ,JI))
-  ENDDO
-ENDDO
-!*     Routing (runoff + drainage)
-!     ----------------------------------------
-!
- CALL ROUTING(ZRUNOFF_TOPD,ZDRAIN_TOPD,KSTEP)
-!
-IF (LHOOK) CALL DR_HOOK('ROUT_DATA_ISBA',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE ROUT_DATA_ISBA
diff --git a/src/SURFEX/routing.F90 b/src/SURFEX/routing.F90
deleted file mode 100644
index 8b8d1874bdf5ba4a33fa045624f19a23b597d3b6..0000000000000000000000000000000000000000
--- a/src/SURFEX/routing.F90
+++ /dev/null
@@ -1,149 +0,0 @@
-!-------------------------------------------------------------------------------
-!     ####################
-      SUBROUTINE ROUTING(PRO,PDR,KSTEP)
-!     ####################
-!
-!!****  *ROUTING*  
-!!
-!!    PURPOSE
-!!    -------
-!     To route the runoff and the exfiltration discharge to the catchment outlet
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* Meteo-France *
-!!      G-M. Saulnier  * LTHE *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   23/11/2005
-!!      M. Le Lay     02/2008 Compatibility with the RESTART option (to update the
-!!                            discharge between two runs)
-!!      Modif B Vincendon 11/2011 : stock managed in thre distinct variables
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_SURF_PAR,        ONLY:XUNDEF
-!
-USE MODD_TOPODYN, ONLY : XTOPD_STEP, NNCAT, XQTOT, NNMC, &
-                         XTIME_TOPD, XQB_RUN, XQB_DR, XTIME_TOPD_DRAIN, NNB_TOPD_STEP
-USE MODD_COUPLING_TOPD, ONLY: XRUN_TOROUT, XDR_TOROUT, NNB_STP_RESTART
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PRO     ! Total water for runoff for each pixel (m3/s)
-!ludo
-REAL, DIMENSION(:,:), INTENT(IN) :: PDR     ! Total water for drainage for each pixel
-INTEGER, INTENT(IN)              :: KSTEP   ! current integration step
-!
-!
-!*      0.2    declarations of local variables
-!
-!
-INTEGER                            :: IUNIT ! unit of discharge files
-INTEGER                            :: JCAT, JJ, JI ! Loop variables
-INTEGER                            :: JSTEP ! current or future integration steps
-REAL, DIMENSION(NNCAT,NNB_TOPD_STEP+NNB_STP_RESTART) :: ZRUN_TOROUT ! stock d'eau dans rivičre pas dans debit ! Kg/m2
-REAL, DIMENSION(NNCAT,NNB_TOPD_STEP+NNB_STP_RESTART) :: ZDR_TOROUT ! stock d'eau dans rivičre pas dans debit ! Kg/m2
- CHARACTER(LEN=3)                   :: YSTEP
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('ROUTING',0,ZHOOK_HANDLE)
-!
-!*       1.0.     Initialization :
-!               --------------
-!
-ZRUN_TOROUT(:,:) = 0.
-ZDR_TOROUT (:,:) = 0.
-!
-DO JCAT=1,NNCAT
-  !
-  !*       2.0    Runoff by geomorpho transfer function
-  !               -------------------------------------
-  DO JJ=1,NNMC(JCAT)
-  !
-    IF ( PRO(JCAT,JJ) > 0.0 .AND. PRO(JCAT,JJ) < XUNDEF ) THEN
-      !
-      JSTEP = INT(XTIME_TOPD(JCAT,JJ) / XTOPD_STEP) + KSTEP 
-      !
-      IF ( JSTEP.LE.NNB_TOPD_STEP ) THEN
-        !
-        XQB_RUN(JCAT,JSTEP) = XQB_RUN(JCAT,JSTEP) + PRO(JCAT,JJ)
-        XQTOT(JCAT,JSTEP) = XQTOT(JCAT,JSTEP) + PRO(JCAT,JJ)
-        !
-      ELSEIF (JSTEP.LE.NNB_TOPD_STEP+NNB_STP_RESTART) THEN
-        !
-        ZRUN_TOROUT(JCAT,JSTEP) = ZRUN_TOROUT(JCAT,JSTEP) + PRO(JCAT,JJ)  !m3
-        !
-      ENDIF
-      !
-    ENDIF
-    !
-  ENDDO
-  !
-  !*       3.0    Drainage by geomorpho transfer function
-  !               -------------------------------------
-  DO JJ=1,NNMC(JCAT)
-    !
-    IF ((PDR(JCAT,JJ) > 0.0).AND.(PDR(JCAT,JJ)<XUNDEF)) THEN
-      !
-      JSTEP = INT(XTIME_TOPD_DRAIN(JCAT,JJ) / XTOPD_STEP) + KSTEP
-      !
-      IF (JSTEP.LE.NNB_TOPD_STEP) THEN
-        !
-        XQB_DR(JCAT,JSTEP) = XQB_DR(JCAT,JSTEP) + PDR(JCAT,JJ)
-        XQTOT(JCAT,JSTEP) = XQTOT(JCAT,JSTEP) + PDR(JCAT,JJ) 
-        !
-      ELSEIF (JSTEP.LE.NNB_TOPD_STEP+NNB_STP_RESTART) THEN
-        !
-        ZRUN_TOROUT(JCAT,JSTEP) = ZRUN_TOROUT(JCAT,JSTEP) + PDR(JCAT,JJ) !m3
-        !
-      ENDIF
-    ENDIF
-    !
-  ENDDO
-  ! 
-  XQB_RUN(JCAT,KSTEP) = XQB_RUN(JCAT,KSTEP) + XRUN_TOROUT(JCAT,KSTEP)
-  XQB_DR(JCAT,KSTEP)  = XQB_DR(JCAT,KSTEP)  + XDR_TOROUT(JCAT,KSTEP)
-  XQTOT(JCAT,KSTEP)   = XQTOT(JCAT,KSTEP) + XRUN_TOROUT(JCAT,KSTEP) + XDR_TOROUT(JCAT,KSTEP)
-  !
-  XRUN_TOROUT(JCAT,:) = XRUN_TOROUT(JCAT,:) + ZRUN_TOROUT(JCAT,:)
-  XDR_TOROUT(JCAT,:)  = XDR_TOROUT(JCAT,:)  + ZDR_TOROUT(JCAT,:)
-  !
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('ROUTING',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE ROUTING
diff --git a/src/SURFEX/sat_area_frac.F90 b/src/SURFEX/sat_area_frac.F90
deleted file mode 100644
index 5b74aad45c2188cf295b1a858cd13df729708573..0000000000000000000000000000000000000000
--- a/src/SURFEX/sat_area_frac.F90
+++ /dev/null
@@ -1,88 +0,0 @@
-!-----------------------------------------------------------------
-!     ########################
-      SUBROUTINE SAT_AREA_FRAC(PDEF,PAS)
-!     ########################
-!
-!!*****    * SAT_AREA_FRAC *
-!
-!!    PURPOSE
-!!    -------
-!    
-!     
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   27/11/2006
-!
-!----------------------------------------------------------------------
-!*       0.      DECLARATIONS
-!                ------------
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, NNMC, XDXT
-USE MODD_COUPLING_TOPD, ONLY : NMASKT
-USE MODD_SURF_ATM_GRID_n, ONLY : XMESH_SIZE
-USE MODD_SURF_PAR,        ONLY : XUNDEF, NUNDEF
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-REAL, DIMENSION(:,:),INTENT(IN)     :: PDEF    ! deficit
-REAL, DIMENSION(:), INTENT(OUT)     :: PAS     !contributive area fraction in Isba meshes
-!
-!*      0.2    declarations of local variables
-INTEGER               :: JJ, JI
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-----------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('SAT_AREA_FRAC',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-! 
-PAS(:)=0.0
-!
-DO JJ=1,NNCAT
-  DO JI=1,NNMC(JJ)
-    IF (PDEF(JJ,JI)==0.0 .AND. NMASKT(JJ,JI)/=NUNDEF) THEN 
-      PAS(NMASKT(JJ,JI)) = PAS(NMASKT(JJ,JI)) + XDXT(JJ)**2
-    ENDIF
-  ENDDO
-ENDDO
-!
-! Calculation of the saturated area ratio in each Isba mesh
-WHERE ((XMESH_SIZE/=0.).AND.(PAS/=XUNDEF))
-  PAS(:) = PAS(:) / XMESH_SIZE(:) 
-ENDWHERE
-!
-IF (LHOOK) CALL DR_HOOK('SAT_AREA_FRAC',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE SAT_AREA_FRAC
diff --git a/src/SURFEX/topd_to_isba.F90 b/src/SURFEX/topd_to_isba.F90
deleted file mode 100644
index de4ba2576fc9e05f7817ddff236edfcf92acb280..0000000000000000000000000000000000000000
--- a/src/SURFEX/topd_to_isba.F90
+++ /dev/null
@@ -1,171 +0,0 @@
-!-----------------------------------------------------------------
-!     ####################
-      SUBROUTINE TOPD_TO_ISBA(KI,KSTEP,GTOPD)
-!     ####################
-!
-!!****  *TOPD_TO_ISBA*  
-!!
-!!    PURPOSE
-!!    -------
-!    
-!     
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* LTHE / Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   09/10/2003
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-!*** SurfEx ***
-USE MODI_UNPACK_SAME_RANK
-!*** Coupling ***
-USE MODI_WRITE_FILE_ISBAMAP
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, NNMC, NNB_TOPD_STEP
-USE MODD_COUPLING_TOPD, ONLY : XWG_FULL, XDTOPT, XWTOPT, XWSUPSAT,&
-                                 NMASKT, XTOTBV_IN_MESH, NNPIX, NFREQ_MAPS_WG
-!
-USE MODD_ISBA_n,          ONLY : XWSAT
-USE MODD_SURF_ATM_GRID_n, ONLY : XMESH_SIZE
-USE MODD_SURF_PAR,        ONLY : XUNDEF
-USE MODD_SURF_ATM_n,      ONLY : NR_NATURE
-USE MODD_ISBA_PAR,        ONLY : XWGMIN
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-!
-INTEGER, INTENT(IN)                 :: KI      ! Grid dimensions
-INTEGER, INTENT(IN)                 :: KSTEP   ! Topodyn current time step
-LOGICAL, DIMENSION(:), INTENT(IN) :: GTOPD     ! 
-!
-!*      0.2    declarations of local variables
-!
-!
-INTEGER                :: JJ, JI          ! loop control 
-INTEGER                :: IUNIT               
-REAL, DIMENSION(KI)    :: ZW              ! TOPODYN water content on ISBA grid (mm)
-REAL, DIMENSION(KI)    :: ZCOUNT          ! TOPO pixel number in an ISBA pixel
-REAL, DIMENSION(KI)    :: ZWSAT_FULL      ! Water content at saturation on the layer 2 
-                                          ! on the full grid
-REAL, DIMENSION(KI)    :: ZWG_OLD
-REAL, DIMENSION(KI)    :: ZDG_FULL
- CHARACTER(LEN=3)       :: YSTEP
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('TOPD_TO_ISBA',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
-ZW(: )= 0.0
-!
-ZWG_OLD(:) = XWG_FULL(:)
-!
-!
-!*       1.     TOPODYN-LAT => ISBA
-!               -------------------
-!*       1.1    mobilizable water
-!               -----------------
-!
-DO JJ=1,NNCAT
-  IF (GTOPD(JJ)) THEN
-    DO JI=1,NNMC(JJ)
-      IF (XDTOPT(JJ,JI) /= XUNDEF) THEN
-        ZW(NMASKT(JJ,JI)) = ZW(NMASKT(JJ,JI)) + XWTOPT(JJ,JI)
-      ENDIF
-    ENDDO
-  ELSE
-    GOTO 10 
-  ENDIF
-ENDDO
-!
-ZCOUNT(:)=REAL(NNPIX(:))
-!
-WHERE (ZCOUNT(:) /= 0.0.AND.XWG_FULL(:)/=XUNDEF)
-  ZW(:) = ZW(:) / ZCOUNT
-ENDWHERE
-!
-!
-! The soil water content is the balanced mean between the soil water content calculated by TOPODYN 
-! and the initial soil water content in each mesh, in function of the area of the mesh occupied by TOPODYN
-WHERE ( XMESH_SIZE(:) - XTOTBV_IN_MESH(:) <= 0.0) ! la maille isba est totalement couverte par des bassins versants
-  XWG_FULL(:) = ZW(:)
-ELSEWHERE (XTOTBV_IN_MESH(:) /= 0.0)
-  XWG_FULL = (XTOTBV_IN_MESH(:)/XMESH_SIZE(:))*ZW(:) + ((XMESH_SIZE(:)-XTOTBV_IN_MESH(:))/XMESH_SIZE(:))*XWG_FULL(:)
-
-ENDWHERE
-!
-XWG_FULL(:) = MAX(XWG_FULL(:),XWGMIN)
-!
-10 CONTINUE 
-!
- CALL UNPACK_SAME_RANK(NR_NATURE,XWSAT(:,2),ZWSAT_FULL)
-!
-IF (.NOT.ALLOCATED(XWSUPSAT)) ALLOCATE(XWSUPSAT(KI))
-XWSUPSAT=0.
-!ludo glace Wsat varie
-WHERE ( XWG_FULL(:) > ZWSAT_FULL(:) .AND. XWG_FULL(:)/=XUNDEF )
-  !ludo calcul sat avant wg
-  XWSUPSAT(:) = XWG_FULL(:) - ZWSAT_FULL(:)
-  XWG_FULL(:) = ZWSAT_FULL(:)
-ENDWHERE
-!
-WHERE (XWG_FULL(:) < XWGMIN .AND. XWG_FULL(:)/=XUNDEF)
-  XWG_FULL(:) = XWGMIN
-ENDWHERE
-!
-IF ( NFREQ_MAPS_WG/=0 .AND. (MOD(KSTEP,NFREQ_MAPS_WG)==0 .OR. KSTEP==NNB_TOPD_STEP) ) THEN
-  ! writing of YSTEP to be able to write maps
-  IF (KSTEP<10) THEN
-    WRITE(YSTEP,'(I1)') KSTEP
-  ELSEIF (KSTEP < 100) THEN
-    WRITE(YSTEP,'(I2)') KSTEP
-  ELSE
-    WRITE(YSTEP,'(I3)') KSTEP
-  ENDIF
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,HFILE='carte_w'//YSTEP,HFORM='FORMATTED',HACTION='WRITE')
-  CALL WRITE_FILE_ISBAMAP(IUNIT,XWG_FULL,KI)
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  !
-ENDIF
-!
-IF (LHOOK) CALL DR_HOOK('TOPD_TO_ISBA',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE TOPD_TO_ISBA
diff --git a/src/SURFEX/topd_to_isba_slope.F90 b/src/SURFEX/topd_to_isba_slope.F90
deleted file mode 100644
index 7771798d368258f2af687e1e2db33c421553dd42..0000000000000000000000000000000000000000
--- a/src/SURFEX/topd_to_isba_slope.F90
+++ /dev/null
@@ -1,108 +0,0 @@
-!-----------------------------------------------------------------
-!     ####################
-      SUBROUTINE TOPD_TO_ISBA_SLOPE(KI)
-!     ####################
-!
-!!****  *TOPD_TO_ISBA*  
-!!
-!!    PURPOSE
-!!    -------
-!    
-!     
-!         
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   12/11/2012
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN,       ONLY : NNCAT, NNMC, XTANB
-USE MODD_COUPLING_TOPD, ONLY : NMASKT,NNPIX
-USE MODD_SURF_ATM_SSO_n, ONLY : XSSO_SLOPE
-USE MODD_SURF_PAR,        ONLY : NUNDEF
-USE MODD_SURF_ATM_n,      ONLY : XNATURE, NSIZE_NATURE, NR_NATURE, &
-                                 NSIZE_FULL, NDIM_NATURE, NDIM_FULL  
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-!
-INTEGER, INTENT(IN)                 :: KI      ! Grid dimensions
-!
-!*      0.2    declarations of local variables
-!
-!
-INTEGER                :: JCAT,JPIX,JJ          ! loop control 
-REAL, DIMENSION(KI)    :: ZCOUNT                ! TOPO pixel number in an ISBA pixel
-                                                ! on the full grid
-REAL, DIMENSION(KI)    :: ZSSO_SLOPE
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('TOPD_TO_ISBA_SLOPE',0,ZHOOK_HANDLE)
-!
-!*        1.0  Compute Mean slope over each ISBA_MESH
-!            ----------------------------------------------------------------------
-!
-!write(*,*) 'pente avt topmodel',MINVAL(XSSO_SLOPE),MAXVAL(XSSO_SLOPE),SUM(XSSO_SLOPE,MASK=XSSO_SLOPE/=XUNDEF)
-!
-ZSSO_SLOPE = XSSO_SLOPE
-!
-ZCOUNT(:) = REAL(NNPIX(:))
-
-WHERE (ZCOUNT /= 0.0)
-   ZSSO_SLOPE = 0.
-ENDWHERE
-!
-DO JCAT=1,NNCAT
-  DO JPIX=1,NNMC(JCAT)
-    IF (NMASKT(JCAT,JPIX) /= NUNDEF) THEN
-      ZSSO_SLOPE(NMASKT(JCAT,JPIX)) = ZSSO_SLOPE(NMASKT(JCAT,JPIX)) + XTANB(JCAT,JPIX)
-    ENDIF
-  ENDDO
-ENDDO
-!
-WHERE (ZCOUNT /= 0.0)
-   ZSSO_SLOPE = ZSSO_SLOPE / ZCOUNT
-ENDWHERE
-!
-XSSO_SLOPE = ZSSO_SLOPE
-!
-!write(*,*) 'pente apres modification',MINVAL(XSSO_SLOPE),MAXVAL(XSSO_SLOPE),COUNT(ZCOUNT/=0.0),SUM(XSSO_SLOPE,MASK=XSSO_SLOPE/=XUNDEF)
-!
-IF (LHOOK) CALL DR_HOOK('TOPD_TO_ISBA_SLOPE',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE TOPD_TO_ISBA_SLOPE
diff --git a/src/SURFEX/topodyn_lat.F90 b/src/SURFEX/topodyn_lat.F90
deleted file mode 100644
index 06d5810b5f84d20c3b46191e6e8e5d9c2280bffc..0000000000000000000000000000000000000000
--- a/src/SURFEX/topodyn_lat.F90
+++ /dev/null
@@ -1,392 +0,0 @@
-!-----------------------------------------------------------------
-!     ###########################
-SUBROUTINE TOPODYN_LAT(PRW,PDEF,PKAPPA,PKAPPAC,GTOPD)
-!,PERROR)
-!     ###########################
-!
-!
-!     PURPOSE
-!     -------
-!     to distribute laterally soil water following topodyn concept
-!
-!
-!     METHOD
-!     ------
-!
-!     EXTERNAL
-!     --------
-!     none
-!
-!
-!     AUTHOR
-!     ------
-!
-!     G.-M. Saulnier * LTHE * 
-!     K. Chancibault * CNRM *
-!
-!     MODIFICATIONS
-!     -------------
-!
-!     Original    12/2003
-!     writing in fortran 90 12/2004
-!------------------------------------------------------------------------------------------
-!
-!*    0.0    DECLARATIONS
-!            ------------
-USE MODD_TOPODYN,       ONLY : NNCAT, NMESHT, XDMAXT, XDXT, XMPARA, NNMC, XCONN, NLINE,&
-                                 XSLOP,  XDAREA
-USE MODD_COUPLING_TOPD, ONLY : XWSTOPT, XDTOPT
-USE MODD_TOPD_PAR,        ONLY : XSTEPK
-!
-USE MODD_SURF_PAR,        ONLY : XUNDEF
-!
-USE MODI_FLOWDOWN
-USE MODI_ABOR1_SFX
-USE MODI_WRITE_FILE_VECMAP
-USE MODI_WRITE_FILE_MAP
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*    0.1   declarations of arguments
-!
-REAL, DIMENSION(:,:), INTENT(IN)   :: PRW
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PDEF
-REAL, DIMENSION(:,:), INTENT(OUT)  :: PKAPPA
-REAL, DIMENSION(:), INTENT(OUT)    :: PKAPPAC
-LOGICAL, DIMENSION(:), INTENT(OUT) :: GTOPD
-!
-!*    0.2   declarations of local variables
-!
-LOGICAL              :: GFOUND  ! logical variable
-REAL                 :: ZSOMME
-REAL                 :: ZM      ! XMPARA in m
-REAL                 :: ZDX     ! XDXT in m
-REAL                 :: ZKVAL, ZKVALMIN, ZKVALMAX
-REAL                 :: ZDAV    ! Averaged deficit (m)
-REAL                 :: ZDAV2   ! Averaged deficit on ZA-ZAS-ZAD (m)
-REAL                 :: ZNDMAXAV,ZNKAV ! temporary averaged maximal deficit and averaged similarity index
-REAL                 :: ZDMAXAV,ZKAV   ! averaged maximal deficit and averaged similarity index
-REAL                 :: ZFUNC
-REAL                 :: ZDIF,ZDIFMIN   ! difference calcul
-REAL                 :: ZNAS, ZNAD     ! temporary saturated and dry relative catchment area
-REAL                 :: ZAS,ZAD        ! saturated and dry relative catchment area 
-REAL                 :: ZA             ! total catchment area
-REAL                 :: ZTMP
-!
-REAL, DIMENSION(NMESHT) :: ZDMAX     ! XDMAXT in m
-REAL, DIMENSION(NMESHT) :: ZRW       ! PRW in m
-REAL, DIMENSION(NMESHT) :: ZDINI     ! initial deficit
-REAL, DIMENSION(NMESHT) :: ZMASK
-REAL, DIMENSION(NMESHT) :: ZKAPPA_PACK, ZDMAX_PACK
-!REAL, DIMENSION(NMESHT) :: ZTMP
-!
-INTEGER              :: J1, J2, JJ !
-INTEGER              :: INKAPPA ! number of steps in similarity index distribution
-INTEGER              :: INPCON  ! number of connected pixels
-INTEGER              :: INAS    ! number of saturated pixels
-INTEGER              :: INAD    ! number of dry pixels
-INTEGER :: I_DIM
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-----------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('TOPODYN_LAT',0,ZHOOK_HANDLE)
-!
-PKAPPA(:,: )= 0.0
-PKAPPAC(:) = 0.
-GTOPD(:) = .TRUE.
-PDEF(:,:) = 0.
-!
-DO JJ = 1,NNCAT
-  !*    0.    Initialisation
-  !           -------------- 
-  ZMASK(:) = 0.0
-  ZRW(:) = 0.0
-  ZDMAX(:) = 0.0
-  INPCON = 0
-  ZDAV = 0.0
-  ZDAV2 = 0.0
-  GFOUND = .FALSE.
-  ZDIFMIN = 99999.
-  !
-  ZRW(:) = PRW(JJ,:) 
-  ZDMAX(:) = XDMAXT(JJ,:)
-  ZDX = XDXT(JJ)
-  ZM = XMPARA(JJ)
-  ZDINI(:) = ZDMAX(:)
-  !
-  !*    0.2   definition of the catchment area concerned by lateral distribution
-  !           ------------------------------------------------------------------
-  !
-  DO J1=1,NNMC(JJ)
-    !
-    IF ( ZRW(J1).GT.0.0 ) THEN
-      ZMASK(J1)=1.0
-    ELSE
-      ZMASK(J1)=0.0
-    ENDIF
-    !
-  ENDDO
-  !
-  CALL FLOWDOWN(NNMC(JJ),ZMASK,XCONN(JJ,:,:),NLINE(JJ,:))
-  !
-  WHERE (ZMASK == 0.0) ZMASK = XUNDEF
-  !
-  !
-  !*    1.    Calcul of hydrological similarity and topographic indexes 
-  !           ---------------------------------------------------------
-  !*    1.1   Calcul of averaged deficit and initialisation of indexes
-  !           --------------------------------------------------------
-  !
-  ZA   = NNMC(JJ) * ZDX**2
-  ZTMP=0.
-  !
-  DO J1=1,NNMC(JJ)
-    !
-    IF (ZMASK(J1)/=XUNDEF) THEN
-      !
-      PKAPPA(JJ,J1) = ZRW(J1) 
-      INPCON = INPCON + 1
-      ZDINI(J1) = ZDMAX(J1) - ZRW(J1)
-      !
-      IF ( ZDINI(J1) <0.0 ) THEN
-        ! WRITE(*,*) J1,'Dini negatif !'
-        ZTMP = ZTMP - ZDINI(J1) !we stock here water above saturation to be
-                                !       distributed among the others pixels
-        ZDINI(J1) = 0.
-      ENDIF
-      !
-      ZDAV = ZDAV + ZDINI(J1)
-      !
-    ELSE
-      !
-      PKAPPA(JJ,J1) = XUNDEF
-      !
-    ENDIF
-    !
-  ENDDO
-  !
-  IF (ZTMP>0.) THEN
-    !write(*,*) COUNT(ZDINI(:)<0.),' pixels avec ZDINI negatif. Volume total :', ZTMP
-    WHERE ( ZDINI(:)>0. ) ZDINI(:) = ZDINI(:)-ZTMP/(COUNT(ZDINI(:)>0.))
-  ENDIF
-  !
-  IF (INPCON >= NNMC(JJ)/1000) THEN
-    !
-    ZDAV = ZDAV / INPCON 
-    ZDAV = ZDAV / ZM
-    !
-    !*    1.2   Propagation of indexes
-    !           ----------------------
-    !
-    CALL FLOWDOWN(NNMC(JJ),PKAPPA(JJ,:),XCONN(JJ,:,:),NLINE(JJ,:))
-    !
-    !*    1.3   Distribution of indexes
-    !           ----------------------
-    !
-    J2=1
-    !
-    DO WHILE ( .NOT.GFOUND .AND. J2.LE.NNMC(JJ) )
-      !
-      IF (ZMASK(J2)/=XUNDEF) THEN
-        !
-        GFOUND = .TRUE.
-        ZKVAL = PKAPPA(JJ,J2) * XDAREA(JJ,J2) / XSLOP(JJ,J2)
-        ZKVAL = LOG(ZKVAL)
-        ZKVALMAX = ZKVAL
-        ZKVALMIN = ZKVAL
-        PKAPPA(JJ,J2) = ZKVAL
-        !
-      ELSE
-        !
-        PKAPPA(JJ,J2) = XUNDEF
-        !
-      ENDIF
-      !
-      J2 = J2 + 1
-      !
-    ENDDO
-    !     
-    DO J1 = J2,NNMC(JJ)
-      !
-      IF (ZMASK(J1)/=XUNDEF) THEN
-        !
-        ZKVAL = PKAPPA(JJ,J1) * XDAREA(JJ,J1) / XSLOP(JJ,J1)
-        ZKVAL = LOG(ZKVAL)
-        !
-        IF (ZKVAL.GT.ZKVALMAX) THEN
-          ZKVALMAX = ZKVAL
-        ELSEIF (ZKVAL.LT.ZKVALMIN) THEN
-          ZKVALMIN = ZKVAL
-        ENDIF
-        !
-        PKAPPA(JJ,J1) = ZKVAL
-        !
-      ELSE
-        !
-        PKAPPA(JJ,J1) = XUNDEF
-        !
-      ENDIF
-      !
-    ENDDO
-    !
-    !*    1.4   Calcul of saturation index
-    !           --------------------------
-    !
-    I_DIM = COUNT( ZMASK(1:NNMC(JJ))/=XUNDEF )
-    ZKAPPA_PACK(:) = XUNDEF
-    ZDMAX_PACK (:) = XUNDEF
-    ZKAPPA_PACK(1:I_DIM) = PACK(PKAPPA(JJ,1:NNMC(JJ)),ZMASK(1:NNMC(JJ))/=XUNDEF)
-    ZDMAX_PACK (1:I_DIM) = PACK(ZDMAX    (1:NNMC(JJ)),ZMASK(1:NNMC(JJ))/=XUNDEF)
-    !
-    INKAPPA = INT((ZKVALMAX - ZKVALMIN) / XSTEPK)
-    !
-    DO J1=1,INKAPPA
-      !
-      ZKVAL = ZKVALMIN + (XSTEPK * (J1-1))
-      INAS = 0
-      INAD = 0
-      ZNDMAXAV = 0.0
-      ZNKAV = 0.0
-      !
-      DO J2=1,I_DIM      
-        !      
-        IF ( ZKAPPA_PACK(J2).GE.ZKVAL ) THEN
-          ! saturated pixel
-          INAS = INAS + 1
-        ELSEIF  (ZKAPPA_PACK(J2).LE.( ZKVAL-(ZDMAX_PACK(J2)/ZM)) ) THEN
-          ! dry pixel
-          INAD = INAD + 1
-          ZNDMAXAV = ZNDMAXAV + ZDMAX_PACK(J2)
-        ELSE
-          ZNKAV = ZNKAV + ZKAPPA_PACK(J2)
-        ENDIF
-        !
-      ENDDO
-      ! 
-      IF (INAD == 0) THEN
-        ZNDMAXAV = 0.0
-      ELSE
-        ZNDMAXAV = ZNDMAXAV /  REAL(INAD)
-      ENDIF
-      !
-      IF ( INPCON == INAS .OR. INPCON == INAD .OR. INPCON == (INAD+INAS)) THEN
-        ZNKAV = 0.0
-      ELSE
-        ZNKAV = ZNKAV / REAL(INPCON - INAD - INAS)
-      ENDIF
-      !
-      IF (INPCON /= 0) THEN
-        ZNAS = REAL(INAS) / REAL(INPCON)
-        ZNAD = REAL(INAD) / REAL(INPCON)
-      ENDIF
-      !
-      ZFUNC = (1 - ZNAS - ZNAD) * ( ZKVAL - ZNKAV )
-      IF (ZM /= 0.) ZFUNC = ZFUNC + (ZNAD * (ZNDMAXAV / ZM))
-      !
-      ZDIF = ABS( ZFUNC - ZDAV )
-      !
-      IF ( ZDIF.LT.ZDIFMIN ) THEN
-        !
-        ZDIFMIN = ZDIF
-        PKAPPAC(JJ) = ZKVAL
-        ZAS = ZNAS
-        ZAD = ZNAD
-        ZDMAXAV = ZNDMAXAV
-        ZKAV = ZNKAV
-        !
-      ENDIF
-      !
-    ENDDO   
-    !
-    !*    2.     Local deficits calculation
-    !            --------------------------
-    !
-    !*    2.1    New averaged deficit on A-Ad-As
-    !            -------------------------------
-    !
-    ZDAV = ZDAV * ZM
-    !
-    IF ( ZAS<1. .AND. ZAD<1. ) THEN
-      !
-      ZDAV2 = (ZDAV - ZDMAXAV * ZAD) / (1 - ZAS - ZAD)
-      !
-    ELSE
-      !
-      IF (ZAS>=1.) WRITE(*,*) 'ALL THE AREA IS SATURATED'
-      IF (ZAD>=1.) WRITE(*,*) 'ALL THE AREA HAS A MAXIMAL DEFICIT'
-      WRITE(*,*) 'ALL THE AREA',ZAS,ZAD
-      CALL ABOR1_SFX("TOPODYN_LAT: ALL THE AREA IS SATURATED OR HAS A MAXIMAL DEFICIT")
-      !
-    ENDIF
-    !
-    !*    2.2    Local deficits 
-    !            --------------
-    !
-    ZSOMME=0.0
-    !ZTMP(:)=XUNDEF
-    !
-    DO J1=1,NNMC(JJ)
-      !
-      IF ( ZMASK(J1)/=XUNDEF ) THEN
-        !
-        IF ( (PKAPPA(JJ,J1).GT.(PKAPPAC(JJ) - ZDMAX(J1)/ZM)) .AND. (PKAPPA(JJ,J1).LT.PKAPPAC(JJ)) ) THEN
-          !
-          PDEF(JJ,J1) = ZM * (ZKAV - PKAPPA(JJ,J1)) + ZDAV2
-          !ZTMP(J1) = 0.5
-          IF (PDEF(JJ,J1) < 0.0) PDEF(JJ,J1) = 0.0
-          !
-        ELSEIF ( PKAPPA(JJ,J1).GE.PKAPPAC(JJ) ) THEN
-          !
-          PDEF(JJ,J1) = 0.0
-          !ZTMP(J1) = 1.
-          !
-        ELSEIF ( PKAPPA(JJ,J1).LE.(PKAPPAC(JJ) - ZDMAX(J1)/ZM) ) THEN
-          !
-          PDEF(JJ,J1) = ZDMAX(J1)
-          !ZTMP(J1) = -1.0
-          !
-        ENDIF
-        !
-        ! nouveau contenu en eau total (m)
-        ZSOMME = ZSOMME + ( XWSTOPT(JJ,J1)*XDTOPT(JJ,J1) - PDEF(JJ,J1) )
-        !
-      ELSE
-        !
-        PDEF(JJ,J1) = ZDMAX(J1)
-        !
-      ENDIF
-      !
-    ENDDO
-    !
-    ! variation du contenu en eau total
-    DO J1=1,NNMC(JJ)
-      !
-      IF (PDEF(JJ,J1)<0.0) THEN
-        WRITE(*,*) 'LAMBDA=',PKAPPA(JJ,J1),'LAMBDAC=',PKAPPAC(JJ)
-      ENDIF
-      !
-    ENDDO
-    !
-    GTOPD(JJ)=.TRUE.
-    !
-  ELSE
-    !
-    !  'Pas de redistribution laterale'
-    GTOPD(JJ)=.FALSE.
-    !
-    PKAPPA(JJ,:) = XUNDEF
-    !  PDEF(JJ,:) = ZDMAX(:) - ZRW(:)
-    PDEF(JJ,:) = ZDINI(:)
-    PKAPPAC(JJ) = XUNDEF
-    !
-  ENDIF
-  ! 
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('TOPODYN_LAT',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE TOPODYN_LAT
-
diff --git a/src/SURFEX/write_budget_coupl_rout.F90 b/src/SURFEX/write_budget_coupl_rout.F90
deleted file mode 100644
index d6f7fc90f3f01ce04300a58af8a5d0b23dfc0d99..0000000000000000000000000000000000000000
--- a/src/SURFEX/write_budget_coupl_rout.F90
+++ /dev/null
@@ -1,137 +0,0 @@
-!
-!     ##########################
-      SUBROUTINE WRITE_BUDGET_COUPL_ROUT
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      L. Bouilloud & B. Vincendon	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original  03/2008 
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-!
-USE MODD_TOPODYN,       ONLY : CCAT, NNCAT, NNB_TOPD_STEP
-USE MODD_BUDGET_COUPL_ROUT
-!
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!      none
-!*      0.2    declarations of local variables
-INTEGER           :: JCAT,JSTP ! loop control
-INTEGER           :: IUNIT     ! file unit numbers
- CHARACTER(LEN=28) :: YFILE     ! file name
- CHARACTER(LEN=40) :: YFORM     ! Writing format
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('WRITE_BUDGET_COUPL_ROUT',0,ZHOOK_HANDLE)
-!
-!*       1.     WRITING BUDGET FILES
-!               ------------
-!
-DO JCAT=1,NNCAT
-  !
-  YFORM='(i6,10f15.1)'
-  YFILE=TRIM('bilan_bv_')//TRIM(CCAT(JCAT))//TRIM('.txt')
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,HFILE=YFILE,HFORM='FORMATTED',HACTION='WRITE')
-  !
-  WRITE(IUNIT,*) '     T','         ',YB_VAR(1),'         ',YB_VAR(2),'         ',&
-                                      YB_VAR(3),'         ',YB_VAR(4),'         ',&
-                                      YB_VAR(5),'         ',YB_VAR(6),'         ',&
-                                      YB_VAR(7),'         ',YB_VAR(8),'         ',&
-                                      YB_VAR(9),'         ',YB_VAR(10)
-  !
-  DO JSTP=1,NNB_TOPD_STEP
-    WRITE(IUNIT,YFORM) JSTP,XB_VAR_BV(JSTP,JCAT,1:10)
-  ENDDO
-  !
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  ! 
-  YFILE=TRIM('bilan_nobv_')//TRIM(CCAT(JCAT))//TRIM('.txt')
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,HFILE=YFILE,HFORM='FORMATTED',HACTION='WRITE')
-  ! 
-  WRITE(IUNIT,*) '     T','         ',YB_VAR(1),'         ',YB_VAR(2),'         ',&
-                                      YB_VAR(3),'         ',YB_VAR(4),'         ',&
-                                      YB_VAR(5),'         ',YB_VAR(6),'         ',&
-                                      YB_VAR(7),'         ',YB_VAR(8),'         ',&
-                                      YB_VAR(9),'         ',YB_VAR(10)
-  !
-  DO JSTP=1,NNB_TOPD_STEP
-    WRITE(IUNIT,YFORM) JSTP,XB_VAR_NOBV(JSTP,JCAT,1:10)
-  ENDDO
-  !
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  !
-  YFORM='(i6,5f15.1)'
-  YFILE=TRIM('bilan_q.txt')
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,HFILE=YFILE,HFORM='FORMATTED',HACTION='WRITE')
-  !
-  WRITE(IUNIT,*) '     T','         ',YB_VARQ(1),'         ',YB_VARQ(2),'         ',&
-                                      YB_VARQ(3),'         ',YB_VARQ(4),'         ',&
-                                      YB_VARQ(5),'         '
-  ! 
-  DO JSTP=1,NNB_TOPD_STEP
-    WRITE(IUNIT,YFORM) JSTP,XB_VAR_Q(JSTP,JCAT,1:5)
-  ENDDO
-  !
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-!  ENDIF 
-  !
-ENDDO !JCAT
-
-YFORM='(i6,10f15.1)'
-YFILE=TRIM('bilan_tot.txt')
-!
- CALL OPEN_FILE('ASCII ',IUNIT,HFILE=YFILE,HFORM='FORMATTED',HACTION='WRITE')
-!   
-WRITE(IUNIT,*) '     T','         ',YB_VAR(1),'         ',YB_VAR(2),'         ',&
-                                    YB_VAR(3),'         ',YB_VAR(4),'         ',&
-                                    YB_VAR(5),'         ',YB_VAR(6),'         ',&
-                                    YB_VAR(7),'         ',YB_VAR(8),'         ',&
-                                    YB_VAR(9),'         ',YB_VAR(10)
-!
-DO JSTP=1,NNB_TOPD_STEP
-  WRITE(IUNIT,YFORM) JSTP,XB_VAR_TOT(JSTP,1:10)
-ENDDO
-!
- CALL CLOSE_FILE('ASCII ',IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('WRITE_BUDGET_COUPL_ROUT',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE WRITE_BUDGET_COUPL_ROUT
diff --git a/src/SURFEX/write_discharge_file.F90 b/src/SURFEX/write_discharge_file.F90
deleted file mode 100644
index af66ffa24e73f6594e7f13cb0b283fff04d80186..0000000000000000000000000000000000000000
--- a/src/SURFEX/write_discharge_file.F90
+++ /dev/null
@@ -1,108 +0,0 @@
-!-----------------------------------------------------------------
-!     #######################
-!
-      SUBROUTINE WRITE_DISCHARGE_FILE(HPROGRAM,HFILE,HFORM,&
-                                      KYEAR,KMONTH,KDAY,KH,KM,PQTOT)
-!     #######################
-!
-!!****  *WRITE_DISCHARGE_FILE*  
-!!
-!!    PURPOSE
-!!    -------
-!     This routine aims at reading topographic files
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!
-!!    
-!!    
-!!
-!!      
-!!    REFERENCE
-!!    ---------
-!!
-!!    
-!!      
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon    * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2006
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN, ONLY  : CCAT, NNCAT, NNB_TOPD_STEP
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
- CHARACTER(LEN=*),      INTENT(IN)  :: HPROGRAM   !
- CHARACTER(LEN=*),      INTENT(IN)  :: HFILE      ! File to be read
- CHARACTER(LEN=*),      INTENT(IN)  :: HFORM      ! Format of the file to be read
-INTEGER, DIMENSION(:), INTENT(IN)  :: KYEAR      ! Year of the beginning of the simulation.
-INTEGER, DIMENSION(:), INTENT(IN)  :: KMONTH     ! Month of the beginning of the simulation.
-INTEGER, DIMENSION(:), INTENT(IN)  :: KDAY       ! Day of the beginning of the simulation.
-INTEGER, DIMENSION(:), INTENT(IN)  :: KH         ! Hour of the beginning of the simulation.
-INTEGER, DIMENSION(:), INTENT(IN)  :: KM         ! Minutes of the beginning of the simulation.
-REAL, DIMENSION(:,:) , INTENT(IN)  :: PQTOT      ! Discharge to be writen
-!
-!
-!*      0.2    declarations of local variables
-!
-INTEGER                   :: JJ,JCAT ! loop control 
-INTEGER                   :: IUNIT       ! Unit of the files
-INTEGER                   :: ILUOUT      ! Unit of the files
-!
- CHARACTER(LEN=28) :: YFILE
- CHARACTER(LEN=40) :: YFORM          ! Writing format
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('WRITE_DISCHARGE_FILE',0,ZHOOK_HANDLE)
-!
-!*       0.3    preparing file openning
-!               ----------------------
-!
- CALL GET_LUOUT(HPROGRAM,ILUOUT)
-YFORM='(I4,A1,I2,A1,I2,A1,I2,A1,I2,A1,F7.2)'
-!
-DO JCAT=1,NNCAT
-  !
-  YFILE = TRIM(CCAT(JCAT))//'_'//TRIM(HFILE)
-  !
-  CALL OPEN_FILE(HPROGRAM,IUNIT,YFILE,HFORM,HACTION='WRITE')
-  !
-  WRITE(IUNIT,*) 'YEAR;MO;DA;HO;MI;',CCAT(JCAT)
-  DO JJ=1,NNB_TOPD_STEP
-    WRITE(IUNIT,YFORM) KYEAR(JJ),';',KMONTH(JJ),';',KDAY(JJ),';',&
-                       KH(JJ)   ,';',KM(JJ)    ,';',PQTOT(JCAT,JJ)
-  ENDDO
-  !
-  CALL CLOSE_FILE(HPROGRAM,IUNIT)
-  !
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('WRITE_DISCHARGE_FILE',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE WRITE_DISCHARGE_FILE
-
diff --git a/src/SURFEX/write_file_isbamap.F90 b/src/SURFEX/write_file_isbamap.F90
deleted file mode 100644
index 7b44b7f734a8b0366831a5d404cf2015dc146f86..0000000000000000000000000000000000000000
--- a/src/SURFEX/write_file_isbamap.F90
+++ /dev/null
@@ -1,101 +0,0 @@
-!-----------------------------------------------------------------
-!     ##########################
-      SUBROUTINE WRITE_FILE_ISBAMAP(KUNIT,PVAR,KI)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   25/01/2005
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-!
-USE MODD_TOPODYN
-USE MODD_SURF_ATM_GRID_n, ONLY : XGRID_PAR
-USE MODD_SURF_PAR,        ONLY : XUNDEF
-!
-USE MODE_GRIDTYPE_CONF_PROJ
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN)             :: KUNIT  ! file unit
-REAL, DIMENSION(:), INTENT(IN)  :: PVAR   ! variable to write in the file
-INTEGER, INTENT(IN)             :: KI    ! Grid dimensions
-!
-!
-!*      0.2    declarations of local variables
-INTEGER                    :: JJ,JI
-INTEGER                    :: JINDEX ! reference number of the pixel
-REAL                       :: ZOUT
-REAL                       :: ZMAX,ZMIN
-REAL, DIMENSION(KI)        :: ZXI, ZYI     ! natural coordinates of ISBA grid (conformal projection)
-REAL, DIMENSION(KI)        :: ZDXI, ZDYI   ! Isba grid resolution in the conformal projection
-INTEGER                    :: IIMAX,IJMAX
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_ISBAMAP',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
- CALL GET_GRIDTYPE_CONF_PROJ(XGRID_PAR,PX=ZXI,PY=ZYI,KIMAX=IIMAX,KJMAX=IJMAX,PDX=ZDXI)
-!
-ZMAX = MAXVAL(PVAR)
-ZMIN = MINVAL(PVAR)
-ZOUT = XUNDEF
-!
-DO JJ=1,5
-  WRITE(KUNIT,*)
-ENDDO
-!
-WRITE(KUNIT,*) ZXI(1)
-WRITE(KUNIT,*) ZYI(1)
-WRITE(KUNIT,*) IIMAX
-WRITE(KUNIT,*) IJMAX
-WRITE(KUNIT,*) ZOUT
-WRITE(KUNIT,*) ZDXI(1)
-WRITE(KUNIT,*) ZMIN
-WRITE(KUNIT,*) ZMAX
-!
-DO JJ = 1,IJMAX
-  DO JI = 1,IIMAX
-    JINDEX = (JJ-1) * IIMAX + JI
-    WRITE(KUNIT,*) PVAR(JINDEX)
-  ENDDO
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_ISBAMAP',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE WRITE_FILE_ISBAMAP
diff --git a/src/SURFEX/write_file_map.F90 b/src/SURFEX/write_file_map.F90
deleted file mode 100644
index a388558860c3a5c1a582e454ecf5f354aacfecda..0000000000000000000000000000000000000000
--- a/src/SURFEX/write_file_map.F90
+++ /dev/null
@@ -1,143 +0,0 @@
-!-----------------------------------------------------------------
-!     ##########################
-      SUBROUTINE WRITE_FILE_MAP(PVAR,HVAR)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   25/01/2005
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN, ONLY : CCAT, NNCAT, NNYC, NNXC, XX0, XY0, XDXT, NLINE, &
-                         XTOPD, XNUL
-!
-USE MODD_SURF_PAR, ONLY : XUNDEF
-!
-USE MODI_GET_LUOUT
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PVAR   ! variable to write in the file
- CHARACTER(LEN=30),    INTENT(IN) :: HVAR   ! end name of the file
-!
-!*      0.2    declarations of local variables
- CHARACTER(LEN=50),DIMENSION(NNCAT) :: CNAME
- CHARACTER(LEN=40)                  :: CFMT
- CHARACTER(*),PARAMETER     :: YPFMT1="('(',I4,'(F10.3,')"
-INTEGER                    :: JWRK1,JJ,JI,JCAT
-INTEGER                    :: IINDEX ! reference number of the pixel
-INTEGER                    :: IUNIT,ILUOUT
-REAL                       :: ZOUT ! pixel not included in the catchment
-REAL                       :: ZMIN,ZMAX
-REAL                       :: ZX1, ZY1, ZX2, ZY2 ! left top and right bottom pixels coordinates
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_MAP',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
- CALL GET_LUOUT('OFFLIN',ILUOUT)
-
-ZOUT = XUNDEF
-!
-DO JCAT=1,NNCAT
-  !
-  CNAME(JCAT) = TRIM(CCAT(JCAT))//TRIM(HVAR)
-  !
-  WRITE(ILUOUT,*) CNAME(JCAT)
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,HFILE=CNAME(JCAT),HFORM='FORMATTED')
-  !
-  !*       1.0    writing header map file
-  !               --------------------------------------
-  !
-  IINDEX = (NNYC(JCAT)-1) * NNXC(JCAT) + 1
-  !
-  ZX1 = XX0(JCAT)
-  ZY1 = XY0(JCAT) + ( (NNYC(JCAT)-1) * XDXT(JCAT) )
-  !
-  ZMIN = MINVAL(PVAR(JCAT,:))
-  ZMAX = MAXVAL(PVAR(JCAT,:),MASK=PVAR(JCAT,:)/=XUNDEF)
-  !
-  DO JJ=1,5
-    WRITE(IUNIT,*)
-  ENDDO
-  !
-  WRITE(IUNIT,*) XX0(JCAT)
-  WRITE(IUNIT,*) XY0(JCAT)
-  WRITE(IUNIT,*) NNXC(JCAT) 
-  WRITE(IUNIT,*) NNYC(JCAT)
-  WRITE(IUNIT,*) ZOUT
-  WRITE(IUNIT,*) XDXT(JCAT)
-  WRITE(IUNIT,*) ZMIN
-  WRITE(IUNIT,*) ZMAX
-  !
-  DO JJ=1,NNYC(JCAT)
-    !
-    DO JI=1,NNXC(JCAT)
-      !
-      IINDEX = (JJ - 1) * NNXC(JCAT) + JI
-      ZX1 = XX0(JCAT) + ((JI-1) * XDXT(JCAT))
-      ZY1 = XY0(JCAT) + ((JJ-1) * XDXT(JCAT))
-      !
-      IF ( XTOPD(JCAT,IINDEX).EQ.XNUL(JCAT) ) THEN
-        !
-        WRITE(IUNIT,*) ZOUT
-        !
-      ELSEIF (NLINE(JCAT,IINDEX)/=0) THEN
-        !
-        WRITE(IUNIT,*) PVAR(JCAT,NLINE(JCAT,IINDEX))
-        !
-      ELSE
-        !
-        WRITE(IUNIT,*) ZOUT
-        !
-      ENDIF
-      !
-    ENDDO
-    !
-  ENDDO
-  !
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  !
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_MAP',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE WRITE_FILE_MAP
diff --git a/src/SURFEX/write_file_masktopd.F90 b/src/SURFEX/write_file_masktopd.F90
deleted file mode 100644
index 6dca1b24f8bb31a8253badb751c0cf086ad4d209..0000000000000000000000000000000000000000
--- a/src/SURFEX/write_file_masktopd.F90
+++ /dev/null
@@ -1,86 +0,0 @@
-!-----------------------------------------------------------------
-!     ##########################
-      SUBROUTINE WRITE_FILE_MASKTOPD(KI)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      B. Vincendon	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   11/2011
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN, ONLY : CCAT, NNCAT
-USE MODD_COUPLING_TOPD, ONLY : NMASKI, NNPIX
-USE MODD_SURF_PAR,        ONLY : XUNDEF
-!
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE MODE_GRIDTYPE_CONF_PROJ
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-INTEGER, INTENT(IN)             :: KI    ! Grid dimensions
-!
-!*      0.2    declarations of local variables
-INTEGER           :: JCAT,JMESH,JPIX
-INTEGER           :: IUNIT
- CHARACTER(LEN=50) :: YNAME
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_MASKTOPD',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
-DO JCAT=1,NNCAT
-  !
-  YNAME = TRIM(CCAT(JCAT))//TRIM('.mask_surf')
-  !
-  CALL OPEN_FILE('ASCII ',IUNIT,YNAME,'FORMATTED',HACTION='WRITE')
-  !
-  DO JMESH=1,KI
-    DO JPIX=1,NNPIX(JMESH)
-      WRITE(IUNIT,*) NMASKI(JMESH,JCAT,JPIX)
-    ENDDO
-  ENDDO
-  !
-  CALL CLOSE_FILE('ASCII ',IUNIT)
-  !
-ENDDO
-!
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_MASKTOPD',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE WRITE_FILE_MASKTOPD
diff --git a/src/SURFEX/write_file_vecmap.F90 b/src/SURFEX/write_file_vecmap.F90
deleted file mode 100644
index a3ff8e379d1d4685ca84a01be999602b654df492..0000000000000000000000000000000000000000
--- a/src/SURFEX/write_file_vecmap.F90
+++ /dev/null
@@ -1,109 +0,0 @@
-!------------------------------------------------------------
-!     ##########################
-      SUBROUTINE WRITE_FILE_VECMAP(PVAR,HVAR,KCAT)
-!     ##########################
-!
-!!
-!!    PURPOSE
-!!    -------
-!        
-!     
-!!**  METHOD
-!!    ------
-!
-!!    EXTERNAL
-!!    --------
-!!
-!!    none
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------ 
-!!      
-!!    REFERENCE
-!!    ---------
-!!     
-!!    AUTHOR
-!!    ------
-!!
-!!      K. Chancibault	* Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!
-!!      Original   25/01/2005
-!-------------------------------------------------------------------------------
-!
-!*       0.     DECLARATIONS
-!               ------------
-!
-USE MODD_TOPODYN
-!
-USE MODD_SURF_PAR, ONLY:XUNDEF
-USE MODI_OPEN_FILE
-USE MODI_CLOSE_FILE
-!
-USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
-USE PARKIND1  ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!*      0.1    declarations of arguments
-!
-REAL, DIMENSION(:),INTENT(IN) :: PVAR   ! variable to write in the file
- CHARACTER(LEN=30), INTENT(IN) :: HVAR   ! end name of the file
-INTEGER,           INTENT(IN) :: KCAT   ! catchment number
-!
-!*      0.2    declarations of local variables
-!
- CHARACTER(LEN=50)          :: CNAME
- CHARACTER(LEN=40)          :: CFMT
-INTEGER                    :: JJ,JI,JK
-INTEGER                    :: IINDEX ! reference number of the pixel
-INTEGER                    :: IUNIT
-REAL                       :: ZOUT ! pixel not included in the catchment
-REAL                       :: ZMIN,ZMAX
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_VECMAP',0,ZHOOK_HANDLE)
-!
-!*       0.     Initialization:
-!               ---------------
-!
-ZOUT = XUNDEF
-ZMIN = MINVAL(PVAR)
-ZMAX = MAXVAL(PVAR)
-!
-CNAME = TRIM(CCAT(KCAT))//TRIM(HVAR)
-!
- CALL OPEN_FILE('ASCII ',IUNIT,HFILE=CNAME,HFORM='FORMATTED')
-!
-DO JI=1,5
-  WRITE(IUNIT,*)
-ENDDO
-!
-WRITE(IUNIT,*) XX0(KCAT)
-WRITE(IUNIT,*) XY0(KCAT)
-WRITE(IUNIT,*) NNXC(KCAT) 
-WRITE(IUNIT,*) NNYC(KCAT)
-WRITE(IUNIT,*) ZOUT
-WRITE(IUNIT,*) XDXT(KCAT)
-WRITE(IUNIT,*) ZMIN
-WRITE(IUNIT,*) ZMAX
-!
-DO JI=1,NNYC(KCAT)
-  DO JK=1,NNXC(KCAT)
-    IINDEX = (JI - 1) * NNXC(KCAT) + JK
-    IF (XTOPD(KCAT,IINDEX).EQ.XNUL(KCAT)) THEN
-      WRITE(IUNIT,*) ZOUT
-    ELSE
-      WRITE(IUNIT,*) PVAR(NLINE(KCAT,IINDEX))
-    ENDIF
-  ENDDO
-ENDDO
-! 
- CALL CLOSE_FILE('ASCII ',IUNIT)
-!
-IF (LHOOK) CALL DR_HOOK('WRITE_FILE_VECMAP',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE WRITE_FILE_VECMAP
-