coupl_topd.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!-----------------------------------------------------------------
!     #####################
      SUBROUTINE coupl_topd (DEC, DC, DMI, PMESH_SIZE, IO, S, K, NK, NP, NPE, &
                             UG, U, 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
!!      09/2013 : Modifications to be able to run with ISBA-DF and more than 1
!!                patch
!!      03/2014: Modif BV : New organisation for first time step (displacement
!!                          from init_coupl_topd)
!!      07/2015: Modif BV : modification of recharge computation
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_diag_n, ONLY : diag_t
USE modd_diag_evap_isba_n, ONLY : diag_evap_isba_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_t, isba_nk_t, isba_np_t, isba_npe_t
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_topd_par, ONLY : nunit
USE modd_topodyn,        ONLY   : nncat, nmesht, nnmc, xmpara, xdmaxt
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, lbudget_topd
                                  !

USE modd_csts,             ONLY : xrholw, xrholi
USE modd_surf_par,         ONLY : xundef, nundef
USE modd_isba_par,         ONLY : xwgmin

USE modd_dummy_exp_profile,ONLY :xf_param, xc_depth_ratio
!
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 modi_dg_dfto3l
USE modi_avg_patch_wg
USE modi_dispatch_wg
USE modi_topd_to_df
USE modi_init_budget_coupl_rout
USE modi_control_water_budget_topd
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(diag_evap_isba_t), INTENT(INOUT) :: DEC
TYPE(diag_t), INTENT(INOUT) :: DC
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMI
REAL, DIMENSION(:), INTENT(IN) :: PMESH_SIZE
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_nk_t), INTENT(INOUT) :: NK
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENt(INOUT) :: NPE
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
!
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(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)           :: ZWM,ZWIM        ! Water content on SurfEx grid after the previous topodyn time step
REAL, DIMENSION(KI)           :: Z_WSTOPI, Z_WFCTOPI
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, ZWSATMAX
LOGICAL, DIMENSION(NNCAT)     :: GTOPD           ! logical variable = true if topodyn_lat runs
INTEGER                       :: JJ, JI, JL, JP    ! loop control 
INTEGER                       :: ILUOUT          ! unit number of listing file
INTEGER                       :: IACT_GROUND_LAYER, IDEPTH, IMASK, ISUM !number of active ground layers
!
REAL, DIMENSION(U%NSIZE_NATURE,3) :: ZWG_3L,ZWGI_3L,ZDG_3L          
 REAL, DIMENSION(U%NSIZE_NATURE)  :: ZMESH_SIZE, ZWSAT         
REAL, DIMENSION(U%NSIZE_NATURE,IO%NGROUND_LAYER,IO%NPATCH) :: ZWG_TMP
REAL, DIMENSION(U%NSIZE_NATURE,IO%NPATCH) :: ZWG, ZDG
REAL, DIMENSION(KI)            :: ZF_PARAM_FULL
REAL, DIMENSION(NNCAT,NMESHT)  :: ZF_PARAMT
!
! Taking several patches into account
INTEGER, DIMENSION(U%NSIZE_NATURE)         :: INB_ACTIVE_PATCH
REAL, DIMENSION(U%NSIZE_NATURE):: ZSUMFRD2, ZSUMFRD3
REAL, DIMENSION(U%NSIZE_NATURE,3) :: ZWG_CTL
 !
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('COUPL_TOPD',0,zhook_handle)
!
CALL GET_LUOUT(HPROGRAM,ILUOUT)
!
!
!*       0.     Initialization:
!               ---------------
zwsatmax=maxval(xwstopi,mask=xwstopi/=xundef)
!
zwg_tmp(:,:,:) = 0.
DO jp = 1,io%NPATCH
  DO jj = 1,np%AL(jp)%NSIZE_P
    imask = np%AL(jp)%NR_P(jj)
    npe%AL(jp)%XWG(jj,:) = max(npe%AL(jp)%XWG(jj,:),xwgmin)
    zwg_tmp(imask,:,jp)  = npe%AL(jp)%XWG(jj,:)
  ENDDO
ENDDO
!
IF (.NOT.ALLOCATED(xwsupsat)) ALLOCATE(xwsupsat(ki))
!
DO jj=1,u%NSIZE_NATURE
 inb_active_patch(jj) = count(s%XPATCH(jj,:)/=0.)
ENDDO
!               ---------------
IF (io%CISBA=='DIF') THEN
 CALL dg_dfto3l(io, np, zdg_3l)
 zwg_3l(:,2)  = dmi%XFRD2_TWG (:)
 zwg_3l(:,3)  = dmi%XFRD3_TWG (:)
 zwgi_3l(:,2) = dmi%XFRD2_TWGI(:)
 zwgi_3l(:,3) = dmi%XFRD3_TWGI(:)

ELSEIF (io%CISBA=='3-L') THEN
 CALL avg_patch_wg(io, np, npe, zwg_3l, zwgi_3l, zdg_3l)
ENDIF
!
zwm(1:ki) = xwg_full(1:ki)
zwim(1:ki) = xwgi_full(1:ki)
!
!*       1.     ISBA => TOPODYN
!               ---------------
!*       1.1    Computation of the useful depth and water for lateral transfers
!               -----------------------------------
!

CALL UNPACK_SAME_RANK(U%NR_NATURE,ZDG_3L(:,2),ZDG2_FULL)
CALL UNPACK_SAME_RANK(U%NR_NATURE,ZDG_3L(:,3),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(U%NR_NATURE,ZWG_3L(:,2),ZWG2_FULL)
CALL UNPACK_SAME_RANK(U%NR_NATURE,ZWG_3L(:,3),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
!
IF (kstep==1) THEN 
  IF (lbudget_topd) CALL init_budget_coupl_rout(dec, dc, dmi, pmesh_size, io, np, npe, u, ki)
 CALL isba_to_topd(xwg_full,xwtopt)
 WHERE (xwtopt == xundef) xwtopt = 0.0
ENDIF
!
!ludo prise en compte glace (pas de glace dans 3e couche)
CALL UNPACK_SAME_RANK(U%NR_NATURE,ZWGI_3L(:,2),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).AND.(xwgi_full/=xundef).AND.(zwim/=xundef))
  zri_wgi = ( (xwgi_full - zwim)  ) * xdtopi!old code
ELSEWHERE
  zri_wgi = 0.0
END WHERE
!CALL UNPACK_SAME_RANK(U%NR_NATURE,DE%XAVG_DWGI(:),ZRI_WGI)
!
WHERE ( xdtopi==xundef ) 
  zri_wgi = 0.0
END WHERE
!
CALL ISBA_TO_TOPD(ZRI_WGI,ZRI_WGIT)
!
!!!!!!!!!!!!!!!!!
!Determination of Wsat, Wfc, Dmax
!!!!!!!!!!!!!!!
!test reservoir top=eau+glace -> pas de modif Wsat et Wfc
  z_wstopi  = xwstopi 
  z_wfctopi = xwfctopi
WHERE ( xwgi_full/=0. .AND.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
!
 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 .AND. xwstopt/=xundef .AND. xdtopt/=xundef)&
                xdmaxfc = (xwstopt - xwfctopt) * xdtopt ! (m)
xdmaxt=xdmaxfc
!WHERE ( XDMAXT >=XUNDEF ) XDMAXT=(MAXVAL(Z_WSTOPI)-MAXVAL(Z_WFCTOPI))*MAXVAL(ZDG2_FULL)

!
!actualisation M
IF( io%CKSAT=='EXP' .OR. io%CKSAT=='SGH' ) THEN
  !ludo test
  xf_param(:) = s%XF_PARAM(:)
  CALL unpack_same_rank(u%NR_NATURE,xf_param(:),zf_param_full)
  CALL isba_to_topd(zf_param_full,zf_paramt)
  !
  !passage de f a M (M=Wsat-Wfc/f)
  !ludo test ksat exp
  WHERE( zf_paramt/=xundef .AND. zf_paramt/=0. ) zf_paramt = (xwstopt-xwfctopt)/zf_paramt
  !
  DO jj=1,nncat
    xmpara(jj) = sum(zf_paramt(jj,:),mask=zf_paramt(jj,:)/=xundef) / nnmc(jj)
  ENDDO
  !
ELSE
  !
  DO jj=1,nncat
    xmpara(jj) = sum( xdmaxfc(jj,:),mask=xdmaxfc(jj,:)/=xundef ) / nnmc(jj) / 4.
  ENDDO
  !
ENDIF
!
!!!!!!!!!!!!!!!
!*       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).AND.(xwg_full/=xundef).AND.(zwm/=xundef))
  zri = ( (xwg_full - zwm)  ) * xdtopi+ zri_wgi
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,GTOPD)
!
CALL PACK_SAME_RANK(U%NR_NATURE,ZAS,XAS_NATURE)
!
!*       3.    Deficit (m) -> water storage (m3/m3) and changing of phase
!               ------------------------------------
!
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
!
!*       4.     TOPODYN => ISBA
!               ---------------
!*       4.1    Calculation of water storage on ISBA grid
!               -----------------------------------------
!
CALL TOPD_TO_ISBA(K, UG, U, KI,KSTEP,GTOPD)!=modif of XWG_FULL from XWTOPT
CALL PACK_SAME_RANK(U%NR_NATURE, (1-XFRAC_D2)*ZWG2_FULL + XFRAC_D2*XWG_FULL, ZWG_3L(:,2))
CALL PACK_SAME_RANK(U%NR_NATURE, (1-XFRAC_D3)*ZWG3_FULL + XFRAC_D3*XWG_FULL, ZWG_3L(:,3))
!
!*       4.2    Budget correction
!  -----------------------------------------
!
 CALL pack_same_rank(u%NR_NATURE,ug%G%XMESH_SIZE,zmesh_size)
 zavg_mesh_size = sum(zmesh_size(:),mask=zmesh_size(:)/=xundef) / count(zmesh_size(:)/=xundef)
!
IF (io%CISBA=='DIF') THEN
 CALL topd_to_df(io, nk, np, npe, zwg_3l)
ELSEIF (io%CISBA=='3-L') THEN
 CALL dispatch_wg(s, np, npe, zwg_3l, zwgi_3l, zdg_3l)
ENDIF
!
DO jp = 1,io%NPATCH
  WHERE(npe%AL(jp)%XWG(:,:)>zwsatmax.AND.npe%AL(jp)%XWG(:,:)/=xundef)
    npe%AL(jp)%XWG(:,:)=zwsatmax
  ENDWHERE
  WHERE(npe%AL(jp)%XWG(:,:)<xwgmin)
    npe%AL(jp)%XWG(:,:)=xwgmin
  ENDWHERE
ENDDO
!
iact_ground_layer=3

IF (io%CISBA=='DIF') THEN

DO jl=2,io%NGROUND_LAYER
  isum = 0
  DO jp = 1,io%NPATCH
    IF (all(npe%AL(jp)%XWG(:,jl)==xundef)) isum = isum + 1
  ENDDO
  IF (isum==io%NPATCH) THEN
         iact_ground_layer=jl-1
         !WRITE(ILUOUT,*) 'IACT_GROUND_LAYER=',IACT_GROUND_LAYER
         EXIT
 ENDIF
ENDDO
!
ENDIF
!
 CALL pack_same_rank(u%NR_NATURE,z_wstopi,zwsat)

zwg(:,:) = 0.
zdg(:,:) = 0.
DO jp = 1,io%NPATCH
  DO jj = 1,np%AL(jp)%NSIZE_P
    imask = np%AL(jp)%NR_P(jj)
    zwg(imask,jp) = npe%AL(jp)%XWG(jj,2)
    zdg(imask,jp) = np%AL(jp)%XDG(jj,2)
  ENDDO
ENDDO
CALL CONTROL_WATER_BUDGET_TOPD(IO, S, U, ZWG_TMP(:,2,:), ZWG, ZDG,&
                               zmesh_size,zavg_mesh_size,zwsat(:))
DO jp = 1,io%NPATCH
  DO jj = 1,np%AL(jp)%NSIZE_P
    imask = np%AL(jp)%NR_P(jj)
    npe%AL(jp)%XWG(jj,2) = zwg(imask,jp)
  ENDDO
ENDDO   
DO jl = 3,iact_ground_layer
  zwg(:,:) = 0.
  zdg(:,:) = 0.
  DO jp = 1,io%NPATCH
    DO jj = 1,np%AL(jp)%NSIZE_P
      imask = np%AL(jp)%NR_P(jj)
      zwg(imask,jp) = npe%AL(jp)%XWG(jj,jl)
      zdg(imask,jp) = np%AL(jp)%XDG(jj,jl)-np%AL(jp)%XDG(jj,jl-1)
    ENDDO
  ENDDO
 CALL control_water_budget_topd(io, s, u, zwg_tmp(:,jl,:), zwg, zdg, &
                                zmesh_size,zavg_mesh_size,zwsat(:))
  DO jp = 1,io%NPATCH
    DO jj = 1,np%AL(jp)%NSIZE_P
      imask = np%AL(jp)%NR_P(jj)
      npe%AL(jp)%XWG(jj,jl) = zwg(imask,jp)
    ENDDO
  ENDDO          

ENDDO
!
DO jp = 1,io%NPATCH
  WHERE(npe%AL(jp)%XWG(:,:)>zwsatmax.AND.npe%AL(jp)%XWG(:,:)/=xundef)
   npe%AL(jp)%XWG(:,:)=zwsatmax
  ENDWHERE
  WHERE(npe%AL(jp)%XWG(:,:)<xwgmin)
   npe%AL(jp)%XWG(:,:)=xwgmin
  ENDWHERE
ENDDO
!
!*      5.0    Total discharge
!              ---------------
!
!*      5.1    Total water for runoff on TOPODYN grid
!              ---------------------------------------
!
!In XAVG_RUNOFFC, the paches have been averaged
CALL UNPACK_SAME_RANK(U%NR_NATURE,DEC%XRUNOFF,ZRUNOFFC_FULL)
CALL UNPACK_SAME_RANK(U%NR_NATURE,XAVG_RUNOFFCM,ZRUNOFFC_FULLM)
!
CALL DIAG_ISBA_TO_ROUT(UG%G%XMESH_SIZE,ZRUNOFFC_FULL,ZRUNOFFC_FULLM,ZRUNOFF_ISBA)
!
xavg_runoffcm(:) = dec%XRUNOFF(:)
!
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
!              ----------------------------------------
!In XAVG_DRAINC, the paches have been average
CALL UNPACK_SAME_RANK(U%NR_NATURE,DEC%XDRAIN*XATOP,ZDRAINC_FULL)
CALL UNPACK_SAME_RANK(U%NR_NATURE,XAVG_DRAINCM*XATOP,ZDRAINC_FULLM)
!
CALL DIAG_ISBA_TO_ROUT(UG%G%XMESH_SIZE,ZDRAINC_FULL,ZDRAINC_FULLM,ZDRAIN_ISBA)
!
WHERE (zdrain_isba==xundef) zdrain_isba=0.
!
xavg_draincm(:)  = dec%XDRAIN(:)
!
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 ',nunit,hfile='carte_surfcont'//hstep,hform='FORMATTED',haction='WRITE')
  CALL write_file_isbamap(ug,nunit,zas,ki)
  CALL close_file('ASCII ',nunit)
ENDIF
!
IF (lhook) CALL dr_hook('COUPL_TOPD',1,zhook_handle)
!
!
END SUBROUTINE coupl_topd