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 (DEC, IO, S, K, NP, NPE, UG, U, HPROGRAM )
!     #######################
!
!!****  *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_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_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_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) :: DEC
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
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=*), INTENT(IN) :: HPROGRAM    ! 
!
!*      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, ZKSAT_NAT       !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(U%NDIM_NATURE,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:
!               ---------------
!
zwg_3l(:,:) = xundef 
zwgi_3l(:,:) = xundef
zdg_3l(:,:) = xundef  
!
ALLOCATE(nmaskt_patch(u%NDIM_NATURE))
!
IF (io%CISBA=='DIF') THEN
 CALL dg_dfto3l(io, np, zdg_3l)
ELSEIF (io%CISBA=='3-L') THEN
 CALL avg_patch_wg(io, np, npe, 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%NDIM_NATURE))
xrunoff_top(:) = dec%XRUNOFF(:)
!
IF(.NOT.ALLOCATED(xavg_runoffcm)) ALLOCATE(xavg_runoffcm(u%NDIM_NATURE))
xavg_runoffcm(:) = dec%XRUNOFF(:)
!
IF(.NOT.ALLOCATED(xavg_draincm )) ALLOCATE(xavg_draincm(u%NDIM_NATURE))
xavg_draincm(:) = dec%XDRAIN(:)
!
!
! Reading masks
 CALL read_file_masktopd(u%NDIM_FULL)
!
!*      2.1     Fraction of SurfEx mesh with TOPMODEL
!               -------------------------------------
!
ALLOCATE(nnbv_in_mesh(u%NDIM_FULL,nncat))
ALLOCATE(xbv_in_mesh(u%NDIM_FULL,nncat))
ALLOCATE(xtotbv_in_mesh(u%NDIM_FULL))
!
xtotbv_in_mesh(:) = 0.0
!
DO jj=1,u%NDIM_FULL
  !
  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%G%XMESH_SIZE(jj)) THEN
    xbv_in_mesh(jj,:) = xbv_in_mesh(jj,:) * ug%G%XMESH_SIZE(jj)/xtotbv_in_mesh(jj)
    xtotbv_in_mesh(jj) = ug%G%XMESH_SIZE(jj)
  ENDIF
ENDDO
!
!*      2.2     Fraction of SurfEx mesh with each catchment
!               -------------------------------------------
!
ALLOCATE(zfrac(u%NDIM_FULL))  ! fraction not covered by catchments
zfrac(:) = ( ug%G%XMESH_SIZE(:)-xtotbv_in_mesh(:) ) / ug%G%XMESH_SIZE(:)
zfrac(:) = min(max(zfrac(:),0.),1.)
!
ALLOCATE(xatop(u%NDIM_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(u%NDIM_FULL))
ALLOCATE(zclay_full(u%NDIM_FULL))
 CALL unpack_same_rank(u%NR_NATURE,k%XSAND(:,2),zsand_full)
 CALL unpack_same_rank(u%NR_NATURE,k%XCLAY(:,2),zclay_full)
!
!ludo prof variable pour tr lat (OK car sol homogene verticalement, faux sinon)
ALLOCATE(zdg2_full(u%NDIM_FULL))
ALLOCATE(zdg3_full(u%NDIM_FULL))
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(u%NDIM_FULL))
zrtop_d2(:) = 0.
!
DO jmesh=1,u%NDIM_FULL
  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(u%NDIM_FULL))
ALLOCATE(xfrac_d3(u%NDIM_FULL))
xfrac_d2(:)=1.
xfrac_d3(:)=0.
!
IF (io%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(u%NDIM_FULL))
WHERE (zdg2_full/=xundef)
  zdg_full = xfrac_d2*zdg2_full + xfrac_d3*(zdg3_full-zdg2_full)
ELSEWHERE
  zdg_full = xundef
END WHERE
!
ALLOCATE(zsandtopi(u%NDIM_FULL))
ALLOCATE(zclaytopi(u%NDIM_FULL))
zsandtopi(:)=0.0
zclaytopi(:)=0.0
ALLOCATE(xdtopi(u%NDIM_FULL))
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(u%NDIM_FULL))
ALLOCATE(xwfctopi(u%NDIM_FULL))
xwstopi(:) = 0.0
xwfctopi(:) = 0.0
xwstopi    = wsat_func_1d(zclaytopi,zsandtopi,io%CPEDOTF)
IF (io%CISBA=='2-L' .OR. io%CISBA=='3-L') THEN
  !  field capacity at hydraulic conductivity = 0.1mm/day
 xwfctopi   = wfc_func_1d(zclaytopi,zsandtopi,io%CPEDOTF)
ELSE IF (io%CISBA=='DIF') THEN
  !  field capacity at water potential = 0.33bar        
 xwfctopi   = w33_func_1d(zclaytopi,zsandtopi,io%CPEDOTF)
END IF
!
!modif ludo test ksat exp
WRITE(iluout,*) 'CKSAT==',io%CKSAT

ALLOCATE(xcstopi(u%NDIM_FULL))
xcstopi(:) = 0.0
IF( io%CKSAT=='SGH' .OR. io%CKSAT=='EXP' ) THEN
  !
  ALLOCATE(zksat_nat(u%NDIM_NATURE))
  zksat_nat(:) = 0.
  ALLOCATE(zksat(u%NDIM_FULL))
  zksat(:) = 0.0
  !
  !valeur patch 1 (idem wsat wfc) a voir cas ou il existe plusieurs patchs 
  CALL unpack_same_rank(np%AL(1)%NR_P,np%AL(1)%XCONDSAT(:,1),zksat_nat)
  CALL unpack_same_rank(u%NR_NATURE,zksat_nat,zksat)
  !
  !passage de definition Ksat(profondeur) en Ksat(deficit)
  WHERE ( zdg_full/=xundef .AND. (xwstopi-xwfctopi/=0.) )
    xcstopi(:) = zksat(:) / (xwstopi(:)-xwfctopi(:))
  END WHERE
  !
  DEALLOCATE(zksat, zksat_nat)
  !
ELSE
  !
  xcstopi(:) = hydcondsat_func(zclaytopi,zsandtopi,io%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(u%NDIM_FULL))
ALLOCATE(xwgi_full(u%NDIM_FULL))
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(u%NDIM_FULL))
ALLOCATE(zwg3_full(u%NDIM_FULL))
!
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(s%XF_PARAM)))
!
! 
!*      5.0      Initial saturated area computation
!               -----------------------------------------------------------
!
ALLOCATE(xas_nature(u%NDIM_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%NR_NATURE, hprogram,u%NDIM_FULL)
!
!*      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