init_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 init_coupl_topd (DGEI, I, UG, U, &
                                  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)
!!      03/2014: Modif BV : New organisation for first time step (displacement
!!                          in coupl_topd)
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
! Modules
!
USE modd_diag_evap_isba_n, ONLY : diag_evap_isba_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_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,nmaskt_patch 
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,  xdmaxt
!
USE modd_surf_par,         ONLY : xundef, nundef

!
! 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 modi_avg_patch_wg
USE modi_dg_dfto3l
!
USE mode_soil
!
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(isba_t), INTENT(INOUT) :: i
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
!
 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   :: 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, DIMENSION(SIZE(I%XWG,1),3)  :: zwg_3l,zwgi_3l,zdg_3l   
!
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:
!               ---------------
ALLOCATE(nmaskt_patch(SIZE(i%XWG,1)))
!
IF (i%CISBA=='DIF') THEN
 CALL dg_dfto3l(i, &
                SIZE(i%XWG,1),zdg_3l)
ELSEIF (i%CISBA=='3-L') THEN
 CALL avg_patch_wg(i, &
                   SIZE(i%XWG,1),zwg_3l,zwgi_3l,zdg_3l)
ENDIF
! 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(u%NSIZE_NATURE))
xrunoff_top(:) = dgei%XAVG_RUNOFFC(:)
!
IF(.NOT.ALLOCATED(xavg_runoffcm)) ALLOCATE(xavg_runoffcm(u%NSIZE_NATURE))
xavg_runoffcm(:) = dgei%XAVG_RUNOFFC(:)
!
IF(.NOT.ALLOCATED(xavg_draincm )) ALLOCATE(xavg_draincm(u%NSIZE_NATURE))
xavg_draincm(:) = dgei%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
  !
  IF (xtotbv_in_mesh(jj)> ug%XMESH_SIZE(jj)) THEN
    xbv_in_mesh(jj,:) = xbv_in_mesh(jj,:) * ug%XMESH_SIZE(jj)/xtotbv_in_mesh(jj)
    xtotbv_in_mesh(jj) = ug%XMESH_SIZE(jj)
  ENDIF
ENDDO
!
!*      2.2     Fraction of SurfEx mesh with each catchment
!               -------------------------------------------
!
ALLOCATE(zfrac(ki))  ! fraction not covered by catchments
zfrac(:) = ( ug%XMESH_SIZE(:)-xtotbv_in_mesh(:) ) / ug%XMESH_SIZE(:)
zfrac(:) = min(max(zfrac(:),0.),1.)
!
ALLOCATE(xatop(u%NSIZE_NATURE)) ! fraction covered by catchments part nature
 CALL pack_same_rank(u%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(u%NR_NATURE,i%XSAND(:,2),zsand_full)
 CALL unpack_same_rank(u%NR_NATURE,i%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(u%NR_NATURE,zdg_3l(:,2),zdg2_full)
 CALL unpack_same_rank(u%NR_NATURE,zdg_3l(:,3),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)*min(xbv_in_mesh(jmesh,jcat)/xtotbv_in_mesh(jmesh),1.)    
    END DO
  ENDIF   
ENDDO
!ZTOP_D2 * D2 < D3 : the depth concerned 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.
!
IF (i%CISBA=='3-L') THEN
 WHERE( zdg2_full/=xundef  ) ! if the depth is < D2
  xfrac_d2(:) = min(1.,zrtop_d2(:))
 END WHERE
 WHERE( zdg2_full/=xundef .AND. zrtop_d2*zdg2_full>zdg2_full  ) ! if the depth is > D2
  xfrac_d3(:) = (zrtop_d2(:)*zdg2_full(:)-zdg2_full(:)) / (zdg3_full(:)-zdg2_full(:))
  xfrac_d3(:) = max(0.,xfrac_d3(:))
 END WHERE
ENDIF
 !
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
xwstopi    = wsat_func_1d(zclaytopi,zsandtopi,i%CPEDOTF)
IF (i%CISBA=='2-L' .OR. i%CISBA=='3-L') THEN
  !  field capacity at hydraulic conductivity = 0.1mm/day
 xwfctopi   = wfc_func_1d(zclaytopi,zsandtopi,i%CPEDOTF)
ELSE IF (i%CISBA=='DIF') THEN
  !  field capacity at water potential = 0.33bar        
 xwfctopi   = w33_func_1d(zclaytopi,zsandtopi,i%CPEDOTF)
END IF
!
!modif ludo test ksat exp
WRITE(iluout,*) 'CKSAT==',i%CKSAT

ALLOCATE(zksat(ki))
zksat(:) = 0.0
ALLOCATE(xcstopi(ki))
xcstopi(:) = 0.0
IF( i%CKSAT=='SGH' .OR. i%CKSAT=='EXP' ) THEN
  !
  !valeur patch 1 (idem wsat wfc) a voir cas ou il existe plusieurs patchs 
  CALL unpack_same_rank(u%NR_NATURE,i%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,i%CPEDOTF)
  !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
  zcoef_anis = 1.
  xcstopt = xcstopt*zcoef_anis
!
!
ALLOCATE(xwg_full(ki))
ALLOCATE(xwgi_full(ki))
ALLOCATE(xwtopt(nncat,nmesht))
xwtopt(:,:) = 0.0
ALLOCATE(xwstopt(nncat,nmesht))
ALLOCATE(xwfctopt(nncat,nmesht))
ALLOCATE(xdmaxfc(nncat,nmesht))
xdmaxfc(:,:) = xundef
ALLOCATE(xdmaxt(nncat,nmesht))
xdmaxt(:,:)=xundef
!
ALLOCATE(zwg2_full(ki))
ALLOCATE(zwg3_full(ki))
!
 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

xwgi_full = 0.
!
!
!*      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 (.NOT.ALLOCATED(xf_param)) ALLOCATE(xf_param(SIZE(i%XF_PARAM)))
!
! 
!*      5.0      Initial saturated area computation
!               -----------------------------------------------------------
!
ALLOCATE(xas_nature(u%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(ug, u, &
                                         hprogram,ki)
!
!*      7.0     deallocate
!               ----------
!
DEALLOCATE(zdg2_full)
DEALLOCATE(zdg3_full)
DEALLOCATE(zwg2_full)
DEALLOCATE(zwg3_full)
!
20 CONTINUE
!
IF (lhook) CALL dr_hook('INIT_COUPL_TOPD',1,zhook_handle)
!
END SUBROUTINE init_coupl_topd