gwf.F90

Go to the documentation of this file.

!     ###################################################################
      SUBROUTINE gwf (TPG,                                                &
                      KLON,KLAT,OPRINT,PTSTEP_RUN,PTSTEP,OMASK,PNUM_AQUI, &
                      PDRAIN,PLEN,PWIDTH,PHC_BED,PTOPO_RIV,PTAUG,         &
                      PAREA,PTRANS,PWEFF,PTABGW_F,PTABGW_H,PHS,           &
                      PHGROUND,PHG_OLD,PSURF_STO,PQGCELL,PWTD,PFWTD,      &
                      PHGHS,PGOUT,PGNEG,                                  &
                      PGSTO_ALL,PGSTO2_ALL,PGIN_ALL,PGOUT_ALL             )
!     ###################################################################
!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!!
!!     
!!
!!    AUTHOR
!!    ------
!!      J.P Vergnes   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original     08/12
!!      09/16   B. Decharme  limit wtd to -1000m
!-------------------------------------------------------------------------------
!
USE modd_trip_grid, ONLY : trip_grid_t
!
USE modd_trip_par
!
USE modi_get_lat_gwf
USE modi_gwf_int
USE modi_gwf_budget
USE modi_gwf_solver
USE modi_gwf_cpl_update
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(trip_grid_t), INTENT(INOUT) :: TPG
!
INTEGER,              INTENT(IN)    :: KLON
INTEGER,              INTENT(IN)    :: KLAT
LOGICAL,              INTENT(IN)    :: OPRINT
!
REAL,                 INTENT(IN)    :: PTSTEP_RUN
REAL,                 INTENT(IN)    :: PTSTEP
!
LOGICAL, DIMENSION(:,:), INTENT(IN) :: OMASK      ! Aquifer mask
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PNUM_AQUI  ! Aquifer ref number
REAL, DIMENSION(:,:), INTENT(IN)    :: PDRAIN     ! Input drainage                [kg/s]
REAL, DIMENSION(:,:), INTENT(IN)    :: PLEN       ! river length                  [m]
REAL, DIMENSION(:,:), INTENT(IN)    :: PWIDTH     ! river widths                  [m]
REAL, DIMENSION(:,:), INTENT(IN)    :: PHC_BED    ! River bed depth               [m]
REAL, DIMENSION(:,:), INTENT(IN)    :: PTOPO_RIV  ! River elevatation             [m]
REAL, DIMENSION(:,:), INTENT(IN)    :: PTAUG      ! ground water transfer time    [s]
REAL, DIMENSION(:,:), INTENT(IN)    :: PAREA      ! Grid-cell area                [m2]
REAL, DIMENSION(:,:), INTENT(IN)    :: PTRANS     ! Transmissivity                [m2/s]
REAL, DIMENSION(:,:), INTENT(IN)    :: PWEFF      ! Effective porosity            [-]
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PTABGW_F   !  Groundwater fraction         [-]
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PTABGW_H   ! Topo height                   [m]
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PHS        ! river height at t             [m]
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PHGROUND   ! water table elevation         [m]
REAL, DIMENSION(:,:), INTENT(INOUT) :: PHG_OLD    ! water table elevation at t-1  [m]
!
REAL, DIMENSION(:,:), INTENT(INOUT), OPTIONAL :: PSURF_STO  ! river channel storage at t    [kg]
!
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL   :: PQGCELL
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL   :: PWTD
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL   :: PFWTD
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL   :: PHGHS
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL   :: PGOUT
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL   :: PGNEG
!
REAL,                 INTENT(OUT), OPTIONAL   :: PGSTO_ALL
REAL,                 INTENT(OUT), OPTIONAL   :: PGSTO2_ALL
REAL,                 INTENT(OUT), OPTIONAL   :: PGIN_ALL
REAL,                 INTENT(OUT), OPTIONAL   :: PGOUT_ALL
!
!*      0.2    declarations of local parameter
!
REAL, PARAMETER       :: ZEPSILON = 1.e-12
INTEGER, PARAMETER    :: IITERMAX = 100
!
REAL, PARAMETER       :: ZGWDZMIN = 10.0    ! Thickness to start to decrease lateral Transmissivity [m]
REAL, PARAMETER       :: ZGWERR   = 0.01    ! Limit of 1cm to limit negatif ground_sto [m]
!
!
!*      0.3    declarations of local variables
!
!
REAL, DIMENSION(KLAT) :: ZLAT
REAL                  :: ZGRID_RES
!
REAL, DIMENSION(KLON,KLAT)          :: ZHDRAIN_RIV      !
REAL, DIMENSION(KLON,KLAT)          :: ZRIVERBED        !
REAL, DIMENSION(KLON,KLAT)          :: ZGWDEEP          !
REAL, DIMENSION(KLON,KLAT)          :: ZCRIV            !
REAL, DIMENSION(KLON,KLAT)          :: ZCC              !
REAL, DIMENSION(KLON,KLAT)          :: ZCR              !
REAL, DIMENSION(KLON,KLAT)          :: ZHCOF            !
REAL, DIMENSION(KLON,KLAT)          :: ZRHS             !
REAL, DIMENSION(KLON,KLAT)          :: ZQDRAIN
REAL, DIMENSION(KLON,KLAT)          :: ZTAUG
REAL, DIMENSION(KLON,KLAT)          :: ZSLOPE
REAL, DIMENSION(KLON,KLAT)          :: ZTRANS
!
REAL                                :: ZEVOL            !
REAL                                :: ZNPTS            ! Number of points in aquifer basins
!
INTEGER                             :: IITER, JLON, JLAT
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
! *     1.  INITIALIZATION
!           --------------
!
IF (lhook) CALL dr_hook('GWF',0,zhook_handle)
!
! * Initialisation variables
!
zhcof(:,:) = xundef
zrhs(:,:) = xundef
zgwdeep(:,:) = 0.0
zcr(:,:) = 0.0
zcc(:,:) = 0.0
zcriv(:,:) = 0.0
zqdrain(:,:) = 0.0
ztrans(:,:) = 0.0
!
! * groundwater mask
!
znpts = REAL(count(omask(:,:)))
!
! * save old hground
!
phg_old(:,:) = phground(:,:)
!
! * Niveau de drainage en riviere
!
zriverbed(:,:) = ptopo_riv(:,:) - phc_bed(:,:)
zhdrain_riv(:,:) = zriverbed(:,:) + min(phc_bed(:,:),phs(:,:))
!
! * grid 
!
 CALL get_lat_gwf(tpg, &
                 klat,zgrid_res,zlat)
!
! * Coefficients nappe/riviere
!
WHERE(omask(:,:))
     zslope(:,:) = min(1.0,max(0.0,phground(:,:)-zriverbed(:,:))/phc_bed(:,:))
     ztaug (:,:) = ptaug(:,:)-(ptaug(:,:)-xday)*zslope(:,:)
     zcriv (:,:) = pwidth(:,:) * plen(:,:)/ztaug(:,:)
ENDWHERE
!
! * Transmissivity
!
WHERE(omask(:,:))
     zslope(:,:) = min(1.0,max(0.0,phground(:,:)-ptopo_riv(:,:)+xgwdzmax-zgwerr)/zgwdzmin)
     ztrans(:,:) = ptrans(:,:)*zslope(:,:)
ENDWHERE
!
 CALL gwf_int(klon,klat,zgrid_res,zlat,omask,pnum_aqui,ztrans,zcr,zcc)
!
! *     2.  ITERATION LOOP
!           --------------
!
zevol    = 1.0
iiter    = 0
!
DO WHILE (zevol>zepsilon.AND.iiter<=iitermax)
!      
  DO jlat=1, klat
     DO jlon=1, klon
!     
        IF(omask(jlon,jlat))THEN
        
!          initialisations des coefficients
           IF(phground(jlon,jlat)<zhdrain_riv(jlon,jlat).AND.phs(jlon,jlat)<=xhsmin)THEN
              zcriv(jlon,jlat) = 0.0
           ENDIF
           IF(phground(jlon,jlat)<zriverbed(jlon,jlat))THEN
              zcriv(jlon,jlat) = 0.0
              zgwdeep(jlon,jlat) = max(0.,(phs(jlon,jlat)-xhsmin))*pwidth(jlon,jlat)*plen(jlon,jlat)/ztaug(jlon,jlat)
           ENDIF
!           
!          formulation des coefficients 
           zhcof(jlon,jlat) = zcriv(jlon,jlat) + (pweff(jlon,jlat)*parea(jlon,jlat)/ptstep_run)
           zrhs(jlon,jlat) = pdrain(jlon,jlat)/xrholw + zgwdeep(jlon,jlat) + zcriv(jlon,jlat)*zhdrain_riv(jlon,jlat)  &
                      + (pweff(jlon,jlat)*parea(jlon,jlat)/ptstep_run)*phg_old(jlon,jlat)
!
        ENDIF
!
     ENDDO
   ENDDO
!
!  approximation
   CALL gwf_solver(klon,klat,znpts,omask,zhcof,zrhs,zcr,zcc,phground,zevol)
!
   iiter = iiter +1
!
ENDDO
!
IF(.NOT.PRESENT(psurf_sto))THEN
  IF (lhook) CALL dr_hook('GWF',1,zhook_handle)
  RETURN
ENDIF
!
! *     3.   WATER BUDGET
!            ------------
!
 CALL gwf_budget(klon,klat,omask,phground,zhdrain_riv, &
                zgwdeep,zcr,zcc,zcriv,pqgcell,zqdrain )
!
DO jlat=1, klat
   DO jlon=1, klon
      IF(omask(jlon,jlat))THEN
         phghs(jlon,jlat) = phground(jlon,jlat)-zhdrain_riv(jlon,jlat)
         pgout(jlon,jlat) = max(0.0,zqdrain(jlon,jlat))
         pgneg(jlon,jlat) = min(0.0,zqdrain(jlon,jlat))
         psurf_sto(jlon,jlat) = psurf_sto(jlon,jlat)+pgneg(jlon,jlat)*ptstep_run
      ELSE
         pgout(jlon,jlat) = pdrain(jlon,jlat)
         pgneg(jlon,jlat) = 0.0
      ENDIF
   ENDDO
ENDDO
!
! *     4.   GLOBAL BUDGET
!            -------------
!
IF(oprint)THEN
!
  pgsto_all  = 0.0
  pgsto2_all = 0.0
  pgin_all   = 0.0
  pgout_all  = 0.0
!
!
  DO jlat=1, klat
     DO jlon=1, klon
        IF(omask(jlon,jlat))THEN
           pgsto_all  = pgsto_all  + pweff(jlon,jlat)*phg_old(jlon,jlat)*xrholw
           pgsto2_all = pgsto2_all + pweff(jlon,jlat)*phground(jlon,jlat)*xrholw
           pgin_all   = pgin_all   + pdrain(jlon,jlat)*ptstep_run/(ptstep*parea(jlon,jlat))
           pgout_all  = pgout_all  + (pgout(jlon,jlat)+pgneg(jlon,jlat)-pqgcell(jlon,jlat)) &
                                   * ptstep_run/(ptstep*parea(jlon,jlat))
        ENDIF
     ENDDO
  ENDDO
!
ENDIF
!
! *     5.   WTD COUPLING
!            ------------
!
pwtd(:,:) = xundef
pfwtd(:,:) = xundef
phg_old(:,:) = xundef
!
 CALL gwf_cpl_update(ptabgw_h,ptabgw_f,omask,ptopo_riv, &
                    phc_bed,phground,phg_old,pwtd,pfwtd)
!
IF (lhook) CALL dr_hook('GWF',1,zhook_handle)
!
END SUBROUTINE gwf