trip_surface_flood.F90

Go to the documentation of this file.

SUBROUTINE trip_surface_flood (KLISTING,PTSTEP,OPRINT,OMASK_FLD,        &
                               PTAB_F,PTAB_H,PTAB_VF,PAREA, PVEL,       &
                               PLEN,PWIDTH,PN_FLOOD,PHC,                &
                               PHS,PSURF_STO,PFLOOD_STO,PSOURCE,        &
                               PFLOOD_STO2,PHFLOOD,PFFLOOD,PFLOOD_LEN,  &
                               PWFLOOD,PQFR,PQRF,PVFIN,PVFOUT,PHSF,     &
                               PFSTO_ALL,PFSTO2_ALL,PSOURCE_ALL,        &
                               PFIN_ALL,PFOUT_ALL,PHF_ALL,PFF_ALL       ) 
!                               
!     ################################################################
!
!!****  *TRIP_SURFACE_FLOOD*  
!!
!!    PURPOSE
!!    -------
!
!     Calculate the flood storage in the next time step based on storages
!     of current time step using the Manning equation.
!     
!     Decharme et al., Clim. Dyn., 2012
!
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!      B. Decharme
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    09/01/14 
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_trip_par,    ONLY : xrholw, xm, xundef
!
USE modi_abort_trip
USE modi_flood_update
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
INTEGER, INTENT(IN)                  :: KLISTING
REAL,    INTENT(IN)                  :: PTSTEP ! Trip timestep value (10800s)
LOGICAL, INTENT(IN)                  :: OPRINT   !Printable budget key
!
LOGICAL, DIMENSION(:,:), INTENT(IN)  :: OMASK_FLD  !Floodplain mask
!
REAL,DIMENSION(:,:,:), INTENT(IN)    :: PTAB_F  ! Flood fraction array
REAL,DIMENSION(:,:,:), INTENT(IN)    :: PTAB_H  ! Topo height array
REAL,DIMENSION(:,:,:), INTENT(IN)    :: PTAB_VF ! Flood volume array
!
REAL,DIMENSION(:,:), INTENT(IN)      :: PAREA      ! Grid-cell area    [m²]
REAL,DIMENSION(:,:), INTENT(IN)      :: PVEL       ! river velocity    [m/s]
REAL,DIMENSION(:,:), INTENT(IN)      :: PLEN       ! river length       [m] 
REAL,DIMENSION(:,:), INTENT(IN)      :: PWIDTH     ! river widths                 [m]
REAL,DIMENSION(:,:), INTENT(IN)      :: PN_FLOOD   ! Manning coeficient over floodplains   [-] (0.1)
REAL,DIMENSION(:,:), INTENT(IN)      :: PHC        ! River bed depth              [m]
REAL,DIMENSION(:,:), INTENT(IN)      :: PHS        ! river channel height [m]
REAL,DIMENSION(:,:), INTENT(IN)      :: PSURF_STO  ! river channel storage at t    [kg]
REAL,DIMENSION(:,:), INTENT(IN)      :: PFLOOD_STO ! Floodplain water storage at t [kg]
REAL,DIMENSION(:,:), INTENT(IN)      :: PSOURCE    ! precip-infiltration-evaporation [kg/s]
!
REAL,DIMENSION(:,:), INTENT(INOUT)   :: PFLOOD_STO2! Floodplain water storage at t+1  [kg]
!
REAL,DIMENSION(:,:), INTENT(OUT)     :: PHFLOOD    ! Floodplain water depth       [m]
REAL,DIMENSION(:,:), INTENT(OUT)     :: PFFLOOD    ! Fraction of flood [-]
REAL,DIMENSION(:,:), INTENT(OUT)     :: PFLOOD_LEN ! Floodplain length along the river [m]
REAL,DIMENSION(:,:), INTENT(OUT)     :: PWFLOOD    ! Floodplain width             [m]
!
REAL,DIMENSION(:,:), INTENT(OUT)     :: PQFR  ! Flood flow to river          [kg/s]
REAL,DIMENSION(:,:), INTENT(OUT)     :: PQRF  ! River flow to floodplain     [kg/s]
REAL,DIMENSION(:,:), INTENT(OUT)     :: PVFIN ! River flow to flood velocity [m/s]
REAL,DIMENSION(:,:), INTENT(OUT)     :: PVFOUT! Flood flow to river velocity [m/s]
REAL,DIMENSION(:,:), INTENT(OUT)     :: PHSF  ! River-Floodplain depth comparison [m] during dt
!
REAL,                 INTENT(OUT)    :: PFSTO_ALL,PFSTO2_ALL,PSOURCE_ALL, &
                                        PFIN_ALL,PFOUT_ALL,PHF_ALL,PFF_ALL 
!                                       Final budget variable
!
!*      0.2    declarations of local variables
!
REAL, PARAMETER                            :: ZLEN_MIN = 1.e-6
!
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZMF
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZSLOPE
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZRADIUS
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZDIST
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZDELTA
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZMF_IN
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZMF_OUT
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZFLD_LEN
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZAREA_SG
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZVRIV
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZVFLD
REAL, DIMENSION(SIZE(PAREA,1),SIZE(PAREA,2)) :: ZVINER
!
REAL    :: ZAREA
!
INTEGER :: ILON, ILAT, JLON, JLAT
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
! * Init
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TRIP_SURFACE_FLOOD',0,zhook_handle)
!
ilon = SIZE(parea,1)
ilat = SIZE(parea,2)
!
pflood_sto2(:,:) = 0.0
phsf(:,:) = 0.0
pvfin(:,:) = 0.0
pvfout(:,:) = 0.0
pqfr(:,:) = 0.0
pqrf(:,:) = 0.0
!
zmf(:,:) = 0.0
zslope(:,:) = 0.0
zradius(:,:) = 0.0
zdist(:,:) = 0.0
zdelta(:,:) = 0.0
zmf_in(:,:) = 0.0
zmf_out(:,:) = 0.0
zfld_len(:,:) = 0.0
!
WHERE(omask_fld(:,:))
   zarea_sg(:,:) = parea(:,:)-(plen(:,:)*pwidth(:,:))
ELSEWHERE
   zarea_sg(:,:) = parea(:,:)
ENDWHERE
!
!-------------------------------------------------------------------------------
! * Update the floodplain storage due to source (Precip inter - LEf - Infil)
!-------------------------------------------------------------------------------
!
WHERE(omask_fld(:,:))
  pflood_sto2(:,:)=pflood_sto(:,:)+ptstep*psource(:,:)
ENDWHERE
!
!------------------------------------------------------------------
! * Update the floodplain geomorphological properties
!------------------------------------------------------------------
! 
 CALL flood_update(ptab_f,ptab_h,ptab_vf,zarea_sg,pflood_sto2, &
                  plen,phflood,pfflood,pflood_len,pwflood     ) 
!   
zfld_len(:,:)=max(zlen_min,pflood_len(:,:))
!
!------------------------------------------------------------------
! * Update the floodplain storage due to inflow or outflow
!------------------------------------------------------------------
!
DO jlat=1,ilat
   DO jlon=1,ilon
!
      IF(omask_fld(jlon,jlat))THEN
!
!        ------------------------------------------------------------------
!        Calculate the potential water mass exchange
!           
         phsf(jlon,jlat)=phs(jlon,jlat)-phc(jlon,jlat)-phflood(jlon,jlat)
!
         zmf(jlon,jlat)=phsf(jlon,jlat)*zfld_len(jlon,jlat)*pwidth(jlon,jlat)*xrholw
         zdist(jlon,jlat)=0.5*(pwidth(jlon,jlat)+pwflood(jlon,jlat))
         zslope(jlon,jlat)=phsf(jlon,jlat)/zdist(jlon,jlat)
!
      ENDIF
!                
!     ------------------------------------------------------------------
!     Floodplain inflow case
!
      IF(zmf(jlon,jlat)>0.0.AND.pfflood(jlon,jlat)<1.0)THEN
!
        zradius(jlon,jlat) = phs(jlon,jlat)-phc(jlon,jlat)
        zradius(jlon,jlat) = exp(xm*log(zradius(jlon,jlat)))
!
        zvfld(jlon,jlat) = zradius(jlon,jlat)*sqrt(zslope(jlon,jlat))/pn_flood(jlon,jlat)
        zvriv(jlon,jlat) = max(pvel(jlon,jlat),zvfld(jlon,jlat))
        zviner(jlon,jlat) = sqrt(zvriv(jlon,jlat)*zvfld(jlon,jlat))
        pvfin(jlon,jlat) = min(zviner(jlon,jlat),zdist(jlon,jlat)/ptstep)        
        pvfout(jlon,jlat) = 0.0
!        
        zmf_in(jlon,jlat) = zmf(jlon,jlat)
        zmf_out(jlon,jlat) = 0.0
!
      ENDIF
!                
!     ------------------------------------------------------------------
!     Floodplain outflow case
!
      IF(zmf(jlon,jlat)<0.0.AND.pfflood(jlon,jlat)>0.0)THEN
!
        zradius(jlon,jlat) = phflood(jlon,jlat)
        zradius(jlon,jlat) = exp(xm*log(zradius(jlon,jlat)))
!
        pvfin(jlon,jlat) = 0.0
        zvfld(jlon,jlat) = zradius(jlon,jlat)*sqrt(-1.0*zslope(jlon,jlat))/pn_flood(jlon,jlat)
        pvfout(jlon,jlat) = min(zvfld(jlon,jlat),zdist(jlon,jlat)/ptstep)
!        
        zmf_in(jlon,jlat) = 0.0
        zmf_out(jlon,jlat) = abs(zmf(jlon,jlat))
!
        zdelta(jlon,jlat) = merge(1.0,0.0,pflood_sto2(jlon,jlat)<zmf_out(jlon,jlat))
!
      ENDIF
!
!     ------------------------------------------------------------------
!     Update the floodplain storage
!
      IF(omask_fld(jlon,jlat).AND.zmf(jlon,jlat)/=0.0)THEN
!
         pqrf(jlon,jlat) = zmf_in(jlon,jlat)*pvfin(jlon,jlat)/zdist(jlon,jlat)
!
         pqfr(jlon,jlat) = (1.0-zdelta(jlon,jlat)) * zmf_out(jlon,jlat)*pvfout(jlon,jlat)/zdist(jlon,jlat) &
                         +      zdelta(jlon,jlat)  * pflood_sto2(jlon,jlat)/ptstep
!
         pflood_sto2(jlon,jlat) = pflood_sto2(jlon,jlat)+ptstep*(pqrf(jlon,jlat)-pqfr(jlon,jlat))      
!
      ENDIF
!
  ENDDO
ENDDO
!
!------------------------------------------------------------------
! * Update the floodplain geomorphological properties
!------------------------------------------------------------------
!
 CALL flood_update(ptab_f,ptab_h,ptab_vf,zarea_sg,pflood_sto2, &
                  plen,phflood,pfflood,pflood_len,pwflood     )
!
!-------------------------------------------------------------------------------
! * Budget calculation
!-------------------------------------------------------------------------------
!
IF(oprint)THEN
!
  pfsto_all   = 0.0
  pfsto2_all  = 0.0
  psource_all = 0.0
  pfin_all    = 0.0
  pfout_all   = 0.0
  phf_all     = 0.0
  pff_all     = 0.0
  zarea       = 0.0
!
  DO jlat=1,ilat
     DO jlon=1,ilon
        IF(omask_fld(jlon,jlat))THEN                
           pfsto_all  = pfsto_all  + pflood_sto(jlon,jlat) / parea(jlon,jlat)
           pfsto2_all = pfsto2_all + pflood_sto2(jlon,jlat) / parea(jlon,jlat)
           pfin_all   = pfin_all   + pqrf(jlon,jlat) / parea(jlon,jlat)
           pfout_all  = pfout_all  + pqfr(jlon,jlat) / parea(jlon,jlat)
           psource_all= psource_all+ psource(jlon,jlat) / parea(jlon,jlat)
           phf_all    = phf_all    + phflood(jlon,jlat) * parea(jlon,jlat)
           pff_all    = pff_all    + pfflood(jlon,jlat) * parea(jlon,jlat)
           zarea      = zarea      + parea(jlon,jlat)
        ENDIF
    ENDDO
  ENDDO
!
  IF(zarea>0.0)THEN
    phf_all = phf_all / zarea
    pff_all = pff_all / zarea
  ENDIF
!
ENDIF
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('TRIP_SURFACE_FLOOD',1,zhook_handle)
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
END SUBROUTINE trip_surface_flood