trip_surface_water.F90

Go to the documentation of this file.

!     #########
SUBROUTINE trip_surface_water (KLISTING,PTSTEP,KGRCN,KSEQ,KNEXTX,KNEXTY,KSEQMAX, &
                               OPRINT,OMASK_VEL,PLEN,PRUNOFF,                    &
                               PVEL,PHS,PSURF_STO,PSURF_STO2,PGOUT,PSIN,PSOUT,   &
                               PAREA,PQFR,PQRF,                                  &
                               PSSTO_ALL,PSSTO2_ALL,PSIN_ALL,PDRUN_ALL,          &
                               PSOUT_ALL,PVEL_ALL,PHS_ALL                        ) 
!     ################################################################
!
!!****  *TRIP_SURFACE_WATER*  
!!
!!    PURPOSE
!!    -------
!
!     Calculate the river storage in the next time step based on the storage of current time step 
!     Where OMASK_VEL=true the Manning equation is used to compute a variable flow velocity.
!
!     
!!**  METHOD
!!    ------
!
!     RK Ordre 4 Rang 4
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!    B. Decharme
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/02/09 
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modn_trip,     ONLY : xcvel
USE modd_trip_par, ONLY : xundef, xm, xvelmin, &
                          xhsmin, xrholw
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abort_trip
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
INTEGER, INTENT(IN)                  :: KLISTING
!
REAL, INTENT(IN)                     :: PTSTEP ! Trip timestep value (10800s)
!
INTEGER, DIMENSION(:,:),INTENT(IN)   :: KGRCN  ! Flow direction (1->8)
INTEGER, DIMENSION(:,:),INTENT(IN)   :: KSEQ   ! River sequence
INTEGER, DIMENSION(:,:),INTENT(IN)   :: KNEXTX ! returns x and y point
INTEGER, DIMENSION(:,:),INTENT(IN)   :: KNEXTY ! of destination grid:
!                                                                    8 1 2
!                                                                    7   3
!                                                                    6 5 4
!
INTEGER, INTENT(IN)                    :: KSEQMAX ! maximum down flow
LOGICAL, INTENT(IN)                    :: OPRINT   !Printable budget key
!
LOGICAL, DIMENSION(:,:), INTENT(IN)    :: OMASK_VEL  ! Variable velocity mask
REAL,    DIMENSION(:,:), INTENT(IN)    :: PLEN       ! river length       [m] 
REAL,    DIMENSION(:,:), INTENT(IN)    :: PAREA      ! Grid-cell area    [m²]
REAL,    DIMENSION(:,:), INTENT(IN)    :: PRUNOFF    ! Surface runoff from ISBA    [kg/s]
REAL,    DIMENSION(:,:), INTENT(IN)    :: PGOUT      ! ground water outflow        [kg/s]
REAL,    DIMENSION(:,:), INTENT(IN)    :: PQFR       ! Flood flow to river         [kg/s]
REAL,    DIMENSION(:,:), INTENT(IN)    :: PQRF       ! River flow to floodplain    [kg/s]
REAL,    DIMENSION(:,:), INTENT(IN)    :: PHS        ! river channel height [m]
REAL,    DIMENSION(:,:), INTENT(IN)    :: PVEL       ! River channel velocity  [m/s]
REAL,    DIMENSION(:,:), INTENT(IN)    :: PSURF_STO  ! river channel storage at t  [kg]
!
REAL,    DIMENSION(:,:), INTENT(INOUT) :: PSURF_STO2 ! river channel storage at t+1[kg]
!
REAL,    DIMENSION(:,:), INTENT(OUT)   :: PSIN  ! Inflow to the surface river reservoir [kg/s]
REAL,    DIMENSION(:,:), INTENT(OUT)   :: PSOUT ! Outflow from the surface river reservoir [kg/s]
!
REAL,                    INTENT(OUT)   :: PSSTO_ALL,PSSTO2_ALL,PSIN_ALL,    &
                                          PDRUN_ALL,PSOUT_ALL,PVEL_ALL,     &
                                          PHS_ALL 
!                                         Final budget variable
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PLEN,1),SIZE(PLEN,2)) :: ZQIN
REAL, DIMENSION(SIZE(PLEN,1),SIZE(PLEN,2)) :: ZRADIUS
REAL, DIMENSION(SIZE(PLEN,1),SIZE(PLEN,2)) :: ZQOUT
REAL, DIMENSION(SIZE(PLEN,1),SIZE(PLEN,2)) :: ZSTOMAX
!
REAL    :: ZAREA
!
INTEGER :: ILON, ILAT, JLON, JLAT, ISEQ, INEXTX, INEXTY
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
! * Init
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TRIP_SURFACE_WATER',0,zhook_handle)
!
ilon = SIZE(plen,1)
ilat = SIZE(plen,2)
!
psurf_sto2(:,:) = 0.0
psin(:,:) = 0.0
psout(:,:) = 0.0
!
zqin(:,:) = 0.0
zradius(:,:) = 0.0
zqout(:,:) = 0.0
zstomax(:,:) = 0.0
!
!-------------------------------------------------------------------------------
! * Sequence loop (optimized computation)
!-------------------------------------------------------------------------------
!
iseq=1
 CALL sequence_loop(iseq)
!
IF(kseqmax>2)THEN
  iseq=2
  CALL sequence_loop(iseq)
ENDIF
!
IF(kseqmax>3)THEN
  DO iseq=3,kseqmax-1
     CALL sequence_loop(iseq)
  ENDDO
ENDIF
!
IF(kseqmax>1)THEN
  iseq=kseqmax
  CALL sequence_loop(iseq)
ENDIF
!
!-------------------------------------------------------------------------------
! * Budget calculation
!-------------------------------------------------------------------------------
!
IF(oprint)THEN
!
  pdrun_all   = 0.0
  pssto_all   = 0.0
  pssto2_all  = 0.0
  psin_all    = 0.0
  psout_all   = 0.0
  pvel_all    = 0.0
  phs_all     = 0.0
  zarea       = 0.0
!
  DO jlat=1,ilat
     DO jlon=1,ilon
        IF(kseq(jlon,jlat)>0)THEN
           pdrun_all  = pdrun_all  + prunoff(jlon,jlat)+pgout(jlon,jlat)+pqfr(jlon,jlat)-pqrf(jlon,jlat)
           pssto_all  = pssto_all  + psurf_sto(jlon,jlat) / parea(jlon,jlat)
           pssto2_all = pssto2_all + psurf_sto2(jlon,jlat) / parea(jlon,jlat)
           psin_all   = psin_all   + zqin(jlon,jlat) / parea(jlon,jlat)
           psout_all  = psout_all  + psout(jlon,jlat) / parea(jlon,jlat)
        ENDIF
        IF(omask_vel(jlon,jlat))THEN
         pvel_all   = pvel_all   + pvel(jlon,jlat) * parea(jlon,jlat)
         phs_all    = phs_all    + phs(jlon,jlat) * parea(jlon,jlat)              
         zarea      = zarea      + parea(jlon,jlat)         
        ENDIF
     ENDDO
  ENDDO
!
  IF(zarea>0.0)THEN
    pvel_all = pvel_all / zarea
    phs_all  = phs_all  / zarea
  ENDIF
!
ENDIF
!
IF (lhook) CALL dr_hook('TRIP_SURFACE_WATER',1,zhook_handle)
!
!-------------------------------------------------------------------------------
 CONTAINS
!-------------------------------------------------------------------------------
!
SUBROUTINE sequence_loop(KNUM)
!
INTEGER, INTENT(IN) :: KNUM
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('TRIP_SURFACE_WATER:SEQUENCE_LOOP',0,zhook_handle)
!
DO jlat=1,ilat
   DO jlon=1,ilon
!      
      IF(kseq(jlon,jlat)==knum)THEN
!
!       ---------------------------------------------------------------------
!       inflow calculation
!
        zqin(jlon,jlat)=zqin(jlon,jlat)+prunoff(jlon,jlat)+pgout(jlon,jlat)+pqfr(jlon,jlat)-pqrf(jlon,jlat)
        psin(jlon,jlat)=zqin(jlon,jlat)
!
!       ------------------------------------------------------------------
!       river channel storage calculation
!
        zstomax(jlon,jlat) = psurf_sto(jlon,jlat)+zqin(jlon,jlat)*ptstep
!
        psurf_sto2(jlon,jlat) = zstomax(jlon,jlat)/(1.0+ptstep*pvel(jlon,jlat)/plen(jlon,jlat))
!
!       -------------------------------------------------------------------
!       supress numerical artifacs
!   
        psurf_sto2(jlon,jlat)=min(zstomax(jlon,jlat),psurf_sto2(jlon,jlat))
!
!       ------------------------------------------------------------------
!       river channel outflow calculation and supress numerical artifacs
!
        zqout(jlon,jlat) = (psurf_sto(jlon,jlat)-psurf_sto2(jlon,jlat))/ptstep+zqin(jlon,jlat)
        psout(jlon,jlat) = max(zqout(jlon,jlat),0.0)
!             
        psurf_sto2(jlon,jlat) = psurf_sto2(jlon,jlat) + (psout(jlon,jlat)-zqout(jlon,jlat))*ptstep        
!
!       ------------------------------------------------------------------
        IF(kgrcn(jlon,jlat)>=1.AND.kgrcn(jlon,jlat)<=8)THEN
          inextx=knextx(jlon,jlat)
          inexty=knexty(jlon,jlat)
          zqin(inextx,inexty)=zqin(inextx,inexty)+psout(jlon,jlat)
        ENDIF
!
      ENDIF
!
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('TRIP_SURFACE_WATER:SEQUENCE_LOOP',1,zhook_handle)
!
END SUBROUTINE sequence_loop
!
!-------------------------------------------------------------------------------
END SUBROUTINE trip_surface_water