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 (DGEI, DGMI, IG, I, 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_evap_isba_n, ONLY : diag_evap_isba_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
USE modd_isba_grid_n, ONLY : isba_grid_t
USE modd_isba_n, ONLY : isba_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) :: dgei
TYPE(diag_misc_isba_t), INTENT(INOUT) :: dgmi
TYPE(isba_grid_t), INTENT(INOUT) :: ig
TYPE(isba_t), INTENT(INOUT) :: i
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, jlayer, jpatch    ! loop control 
INTEGER                       :: iluout          ! unit number of listing file
INTEGER                       :: iact_ground_layer, idepth !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,I%NGROUND_LAYER,I%NPATCH) :: zwg_tmp
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)
i%XWG(:,:,:)     = max(i%XWG(:,:,:),xwgmin)
zwg_tmp(:,:,:) = i%XWG(:,:,:)
IF (.NOT.ALLOCATED(xwsupsat)) ALLOCATE(xwsupsat(ki))
!
DO jj=1,u%NSIZE_NATURE
 inb_active_patch(jj)=count(i%XPATCH(jj,:)/=0.)
ENDDO
!               ---------------
IF (i%CISBA=='DIF') THEN
 CALL dg_dfto3l(i, &
                SIZE(i%XWG,1),zdg_3l)
 zwg_3l(:,2)=dgmi%XFRD2_TWG(:)
 zwg_3l(:,3)=dgmi%XFRD3_TWG(:)
 zwgi_3l(:,2)=dgmi%XFRD2_TWGI(:)
 zwgi_3l(:,3)=dgmi%XFRD3_TWGI(:)

ELSEIF (i%CISBA=='3-L') THEN
 CALL avg_patch_wg(i, &
                   SIZE(i%XWG,1),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(dgei, dgmi, ig, i, 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,DGEI%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( i%CKSAT=='EXP' .OR. i%CKSAT=='SGH' ) THEN
  !ludo test
  xf_param(:) = i%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(i, 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%XMESH_SIZE,zmesh_size)
 zavg_mesh_size = sum(zmesh_size(:),mask=zmesh_size(:)/=xundef) / count(zmesh_size(:)/=xundef)
!
IF (i%CISBA=='DIF') THEN
 CALL topd_to_df(i, &
                 SIZE(zwg_3l,1),zwg_3l)
ELSEIF (i%CISBA=='3-L') THEN
 CALL dispatch_wg(i, &
                  SIZE(zwg_3l,1),zwg_3l,zwgi_3l,zdg_3l)
ENDIF
!
WHERE(i%XWG(:,:,:)>zwsatmax.AND.i%XWG(:,:,:)/=xundef)
 i%XWG(:,:,:)=zwsatmax
ENDWHERE
WHERE(i%XWG(:,:,:)<xwgmin)
 i%XWG(:,:,:)=xwgmin
ENDWHERE
!
iact_ground_layer=3

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

DO jlayer=2,i%NGROUND_LAYER
 IF (all(i%XWG(:,jlayer,:)==xundef))THEN
         iact_ground_layer=jlayer-1
         !WRITE(ILUOUT,*) 'IACT_GROUND_LAYER=',IACT_GROUND_LAYER
         EXIT
 ENDIF
ENDDO
!
ENDIF
!
 CALL pack_same_rank(u%NR_NATURE,z_wstopi,zwsat)
 CALL control_water_budget_topd(i, u, &
                               zwg_tmp(:,2,:), i%XWG(:,2,:),i%XDG(:,2,:),&
           zmesh_size,zavg_mesh_size,zwsat(:))
DO jlayer = 3,iact_ground_layer
 CALL control_water_budget_topd(i, u, &
                               zwg_tmp(:,jlayer,:), i%XWG(:,jlayer,:),&
              i%XDG(:,jlayer,:)-i%XDG(:,jlayer-1,:),&
              zmesh_size,zavg_mesh_size,zwsat(:))
ENDDO
!
WHERE(i%XWG(:,:,:)>zwsatmax.AND.i%XWG(:,:,:)/=xundef)
 i%XWG(:,:,:)=zwsatmax
ENDWHERE
WHERE(i%XWG(:,:,:)<xwgmin)
 i%XWG(:,:,:)=xwgmin
ENDWHERE
!
!*      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,dgei%XAVG_RUNOFFC,zrunoffc_full)
 CALL unpack_same_rank(u%NR_NATURE,xavg_runoffcm,zrunoffc_fullm)
!
 CALL diag_isba_to_rout(ug, &
                       zrunoffc_full,zrunoffc_fullm,zrunoff_isba)
!
xavg_runoffcm(:) = dgei%XAVG_RUNOFFC(:)
!
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,dgei%XAVG_DRAINC*xatop,zdrainc_full)
 CALL unpack_same_rank(u%NR_NATURE,xavg_draincm*xatop,zdrainc_fullm)
!
 CALL diag_isba_to_rout(ug, &
                       zdrainc_full,zdrainc_fullm,zdrain_isba)
!
WHERE (zdrain_isba==xundef) zdrain_isba=0.
!
xavg_draincm(:)  = dgei%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 ',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