recharge_surf_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 recharge_surf_topd(PHI,PHT,KI)
!     #######################
!
!!****  *RECHARGE_SURF_TOPD*  
!!
!!    PURPOSE
!!    -------
!
!     
!         
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!
!!      K. Chancibault  * LTHE / Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!
!!      Original   12/2003
!!                 03/2014 (B. Vincendon) use of the number of pixels included in a mesh and a watershed
!! 
!!    WARNING
!!    ----------------
!!     WFC is the threshold for deficits 
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
USE modd_coupling_topd, ONLY: nmaski, nmaskt, xwfctopt, xdmaxfc, xwtopt,&
                                xdtopt, xwstopt, nnpix, nnbv_in_mesh
USE modd_topodyn,       ONLY: nncat, xdmaxt
!
USE modd_surf_par,        ONLY: nundef,xundef
!
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
INTEGER, INTENT(IN) :: KI    ! Grid dimensions
REAL, DIMENSION(:), INTENT(INOUT)   :: PHI   ! water content variation since last time step from ISBA (m)
REAL, DIMENSION(:,:), INTENT(OUT)   :: PHT   ! water content variation to provide to TOPODYN to be distributed (m) 
!
!*      0.2    declarations of local variables
!
!
LOGICAL, DIMENSION(NNCAT,SIZE(NMASKI,3)) :: GTEST
INTEGER            :: J1,J2,J3,J4     ! loop control 
INTEGER            :: INBSAT, INBALL
!
REAL                           :: ZREST            ! m
REAL                           :: ZWNEW            ! m3/m3
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('RECHARGE_SURF_TOPD',0,zhook_handle)
!
!*       0.     Initialization:
!               ---------------
!
!*       1.     ISBA => TOPODYN-LAT
!               -------------------
!
pht(:,:)=0.
!
DO j3 = 1,ki
  !
  !The water content is lower than the previous one : this case is dealed with first to fasten the computation time.
  IF (phi(j3) <= 0.0) THEN
    !
    DO j1 = 1,nncat
      !
      j4 = 1
      j2 = nmaski(j3,j1,j4)
      DO WHILE (j2 /= nundef .AND. j4<=nnbv_in_mesh(j3,j1) )
        !
        IF ( nmaskt(j1,j2) /= nundef ) THEN
          !
          zwnew = xwtopt(j1,j2) + phi(j3) / xdtopt(j1,j2)
          !
          IF ( zwnew >= xwfctopt(j1,j2) ) THEN
            !
            ! on reste au-dessus de la capacite au champ, malgre l'assechement
            IF (xdmaxfc(j1,j2)/=xundef) xdmaxt(j1,j2) = xdmaxfc(j1,j2) 
            pht(j1,j2) = (zwnew - xwfctopt(j1,j2)) * xdtopt(j1,j2)
            !
          ELSE ! on passe au-dessous de la capacite au champ
            !
            IF (xwstopt(j1,j2)/=xundef) &
                    xdmaxt(j1,j2) = (xwstopt(j1,j2) - zwnew) * xdtopt(j1,j2)
            pht(j1,j2) = 0.0
            !
          ENDIF
          !
          j4 = j4+1
          IF ( j4<=SIZE(nmaski,3) ) j2 = nmaski(j3,j1,j4)
          !
        ELSE ! pixel non défini dans Isba
          !
          xdmaxt(j1,j2) = 0.0
          pht = 0.0
          !
        ENDIF
        !
      ENDDO
      !
    ENDDO
    !
  ELSE ! recharge > 0.0
    !
    zrest=1.
    gtest(:,:)=.true.
    !
    DO WHILE ( zrest>0.0 )
      !
      zrest=0.0
      !
      DO j1=1,nncat
        !
        j4=1
        j2=nmaski(j3,j1,j4)
        !
        DO WHILE ( j2/=nundef .AND. j4<=nnbv_in_mesh(j3,j1) )
          !
          IF ( gtest(j1,j4) .AND. nmaskt(j1,j2)/=nundef ) THEN
            !
            zwnew = xwtopt(j1,j2) + phi(j3) / xdtopt(j1,j2)
            !
            IF ( xwtopt(j1,j2) == xwstopt(j1,j2) ) THEN ! pixel déjà saturé              !
              XDMAXT(J1,J2) = 0.0
              PHT(J1,J2) = 0.0
              ZREST = ZREST + PHI(J3)
              GTEST(J1,J4) = .FALSE.
              !
            ELSE IF ( ( XWSTOPT(J1,J2) - XWTOPT(J1,J2) ) * XDTOPT(J1,J2) <= PHI(J3) ) THEN
              !
              ! pixel va se saturer
              XDMAXT(J1,J2) = XDMAXFC(J1,J2)
              PHT(J1,J2) = ( XWSTOPT(J1,J2) - XWFCTOPT(J1,J2) ) * XDTOPT(J1,J2)
              ZREST = ZREST + PHI(J3) - PHT(J1,J2)
              GTEST(J1,J4)=.FALSE.
              !
            ELSE IF ( XWTOPT(J1,J2) < XWFCTOPT(J1,J2) ) THEN 
              !
              ! en dessous de la capacité au champ avant d'ajouter la recharge
              IF ( (XWTOPT(J1,J2) + PHI(J3)/XDTOPT(J1,J2)) <= XWFCTOPT(J1,J2) ) THEN 
                !
                ! en dessous de la capacité au champ avec la recharge 
                XDMAXT(J1,J2) = ( XWSTOPT(J1,J2) - ZWNEW ) * XDTOPT(J1,J2)
                PHT(J1,J2) = 0.0
                !
              ELSE ! au-dessus de la capacité au champ avec la recharge
                !
                XDMAXT(J1,J2) = XDMAXFC(J1,J2)
                PHT(J1,J2) = ( ZWNEW - XWFCTOPT(J1,J2) ) * XDTOPT(J1,J2)
                !
              ENDIF
              !
            ELSE ! au-dessus de la capacité au champ avant d'ajouter la recharge
              !
              XDMAXT(J1,J2) = XDMAXFC(J1,J2)
              PHT(J1,J2) = ( ZWNEW - XWFCTOPT(J1,J2) ) * XDTOPT(J1,J2)
              !
            ENDIF
            !
          ELSE IF ( NMASKT(J1,J2)==NUNDEF ) THEN! pixel non défini dans Isba
            !
            XDMAXT(J1,J2) = 0.0
            PHT = 0.0
            !
          ENDIF
          !
          J4 = J4+1
          IF ( J4<=SIZE(NMASKI,3) ) J2 = NMASKI(J3,J1,J4)
          !
        ENDDO
        !
      ENDDO
      !
      IF ( ZREST/=0.0 ) THEN
        !
        INBSAT=COUNT(.NOT.GTEST) !nb de pixels saturés avec ou sans la recharge
        !
        IF ( INBSAT == NNPIX(J3) ) THEN
          !
          IF (NNPIX(J3) > 400 ) THEN
            WRITE(*,*) 'MAILLE NUM=',J3, 'nb pix tot=',NNPIX(J3)
           ! CALL ABOR1_SFX("RECHARGE_SURF_TOPD: TOO MANY PIXELS SATURATED ")
          ELSE
            ZREST=0.0
          ENDIF
          !
        ELSE
          !
          PHI(J3) = PHI(J3) + ( ZREST / (NNPIX(J3) - INBSAT) ) ! nouvelle recharge à distribuer
          !
        ENDIF
      ENDIF
      !
    ENDDO
    !
  ENDIF
  !
ENDDO
!
IF (LHOOK) CALL DR_HOOK('RECHARGE_SURF_TOPD',1,ZHOOK_HANDLE)
!
END SUBROUTINE RECHARGE_SURF_TOPD

              !
              xdmaxt(j1,j2) = 0.0
              pht(j1,j2) = 0.0
              zrest = zrest + phi(j3)
              gtest(j1,j4) = .false.
              !
            ELSE IF ( ( xwstopt(j1,j2) - xwtopt(j1,j2) ) * xdtopt(j1,j2) <= phi(j3) ) THEN
              !
              ! pixel va se saturer
              xdmaxt(j1,j2) = xdmaxfc(j1,j2)
              pht(j1,j2) = ( xwstopt(j1,j2) - xwfctopt(j1,j2) ) * xdtopt(j1,j2)
              zrest = zrest + phi(j3) - pht(j1,j2)
              gtest(j1,j4)=.false.
              !
            ELSE IF ( xwtopt(j1,j2) < xwfctopt(j1,j2) ) THEN 
              !
              ! en dessous de la capacité au champ avant d'ajouter la recharge
              IF ( (xwtopt(j1,j2) + phi(j3)/xdtopt(j1,j2)) <= xwfctopt(j1,j2) ) THEN 
                !
                ! en dessous de la capacité au champ avec la recharge 
                xdmaxt(j1,j2) = ( xwstopt(j1,j2) - zwnew ) * xdtopt(j1,j2)
                pht(j1,j2) = 0.0
                !
              ELSE ! au-dessus de la capacité au champ avec la recharge
                !
                xdmaxt(j1,j2) = xdmaxfc(j1,j2)
                pht(j1,j2) = ( zwnew - xwfctopt(j1,j2) ) * xdtopt(j1,j2)
                !
              ENDIF
              !
            ELSE ! au-dessus de la capacité au champ avant d'ajouter la recharge
              !
              xdmaxt(j1,j2) = xdmaxfc(j1,j2)
              pht(j1,j2) = ( zwnew - xwfctopt(j1,j2) ) * xdtopt(j1,j2)
              !
            ENDIF
            !
          ELSE IF ( nmaskt(j1,j2)==nundef ) then! pixel non défini dans Isba
            !
            xdmaxt(j1,j2) = 0.0
            pht = 0.0
            !
          ENDIF
          !
          j4 = j4+1
          IF ( j4<=SIZE(nmaski,3) ) j2 = nmaski(j3,j1,j4)
          !
        ENDDO
        !
      ENDDO
      !
      IF ( zrest/=0.0 ) THEN
        !
        inbsat=count(.NOT.gtest) !nb de pixels saturés avec ou sans la recharge
        !
        IF ( inbsat == nnpix(j3) ) THEN
          !
          IF (nnpix(j3) > 400 ) THEN
            WRITE(*,*) 'MAILLE NUM=',j3, 'nb pix tot=',nnpix(j3)
           ! CALL ABOR1_SFX("RECHARGE_SURF_TOPD: TOO MANY PIXELS SATURATED ")
          ELSE
            zrest=0.0
          ENDIF
          !
        ELSE
          !
          phi(j3) = phi(j3) + ( zrest / (nnpix(j3) - inbsat) ) ! nouvelle recharge à distribuer
          !
        ENDIF
      ENDIF
      !
    ENDDO
    !
  ENDIF
  !
ENDDO
!
IF (lhook) CALL dr_hook('RECHARGE_SURF_TOPD',1,zhook_handle)
!
END SUBROUTINE recharge_surf_topd