topodyn_lat.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 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, xlambda
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.
zas=0.
zad=1.
!
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.
  !
  zkav = 0.
  !
  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)>0.0 .AND. zrw(j1)/=xundef) 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) * exp(xlambda(jj,j2))
        !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) * exp(xlambda(jj,j1))
!       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. .AND. (zas + 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