gwf_cpl_update.F90

Go to the documentation of this file.

      SUBROUTINE gwf_cpl_update(PTABGW_H,PTABGW_F,OMASK_GW,PTOPO_RIV, &
                                PHC_BED,PHGROUND,PHG_OLD,PWTD,PFWTD   )
!     ##########################################################################
!
!!****  *GWF_CPL_UPDATE*  
!!
!!    PURPOSE
!!    -------
!
!     update groundwater diag
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!      B. Decharme
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/11/06 
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modn_trip, ONLY : lgwsubf, xgwsubd
!
USE modi_abort_trip
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL,    DIMENSION(:,:,:), INTENT(IN)    :: PTABGW_H
REAL,    DIMENSION(:,:,:), INTENT(IN)    :: PTABGW_F
LOGICAL, DIMENSION(:,:),   INTENT(IN)    :: OMASK_GW
REAL,    DIMENSION(:,:),   INTENT(IN)    :: PTOPO_RIV
REAL,    DIMENSION(:,:),   INTENT(IN)    :: PHC_BED
REAL,    DIMENSION(:,:),   INTENT(IN)    :: PHGROUND
!
REAL,    DIMENSION(:,:),   INTENT(INOUT) :: PHG_OLD
REAL,    DIMENSION(:,:),   INTENT(OUT)   :: PWTD
REAL,    DIMENSION(:,:),   INTENT(OUT)   :: PFWTD
!
!*      0.2    declarations of local variables
!
INTEGER, DIMENSION(SIZE(PTABGW_H,1),SIZE(PTABGW_H,2)) :: ISUP
INTEGER, DIMENSION(SIZE(PTABGW_H,1),SIZE(PTABGW_H,2)) :: IINF
LOGICAL, DIMENSION(SIZE(PTABGW_H,1),SIZE(PTABGW_H,2)) :: LMASK
REAL,    DIMENSION(SIZE(PTABGW_H,1),SIZE(PTABGW_H,2)) :: ZSLOPE
REAL,    DIMENSION(SIZE(PTABGW_H,1),SIZE(PTABGW_H,2)) :: ZHGROUND
!
INTEGER         :: ILON, ILAT, JLON, JLAT, JFRAC, INFRAC
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('GWF_CPL_UPDATE',0,zhook_handle)
!
!-------------------------------------------------------------------------------
! * Init
!-------------------------------------------------------------------------------
!
ilon   = SIZE(ptabgw_h,1)
ilat   = SIZE(ptabgw_h,2)
infrac = SIZE(ptabgw_h,3)
!
zhground(:,:) = 0.0
zslope(:,:) = 0.0

lmask(:,:) = (omask_gw(:,:).AND.phground(:,:)/=phg_old(:,:))
!
!-------------------------------------------------------------------------------
! * Evolution of water table depth
!-------------------------------------------------------------------------------
!
WHERE(lmask(:,:))
      pwtd(:,:) = phground(:,:)-ptopo_riv(:,:)
ENDWHERE
!
!-------------------------------------------------------------------------------
! * Evolution of the fraction of water table to rise
!-------------------------------------------------------------------------------
!
IF(.NOT.lgwsubf)THEN
!
! No sub-grid fraction
!
  WHERE(omask_gw(:,:))
        pfwtd(:,:) = 1.0
  ENDWHERE
!
ELSE
!
! Adjust depth before computed sub-grid fraction
!
  WHERE(lmask(:,:))
    zslope(:,:) = min(1.0,max(0.0,phground(:,:)-(ptopo_riv(:,:)-phc_bed(:,:)))/phc_bed(:,:))
    zhground(:,:) = phground(:,:) + zslope(:,:)*xgwsubd
  ENDWHERE
!
! Compute sub-grid fraction as in Verges et al., JGR, 2014
!
  WHERE(lmask(:,:).AND.zhground(:,:)<=ptabgw_h(:,:,1))
        pfwtd(:,:) = min(1.0,ptabgw_f(:,:,1))
        lmask(:,:) = .false.
  ELSEWHERE(lmask(:,:).AND.zhground(:,:)>=ptabgw_h(:,:,infrac))
        pfwtd(:,:) = min(1.0,ptabgw_f(:,:,infrac))
        lmask(:,:) = .false.
  ENDWHERE
!
  isup(:,:)=0
  iinf(:,:)=0
!
  DO jlat=1,ilat
     DO jlon=1,ilon 
        IF(lmask(jlon,jlat))THEN
          DO jfrac=1,infrac-1
             IF(zhground(jlon,jlat)>=ptabgw_h(jlon,jlat,jfrac))THEN
               isup(jlon,jlat)=jfrac+1
               iinf(jlon,jlat)=jfrac
             ENDIF          
          ENDDO
          IF(iinf(jlon,jlat)==0.or.isup(jlon,jlat)==0)then
            WRITE(6,*)'IINF,ISUP,JLON,JLAT',iinf(jlon,jlat),isup(jlon,jlat),jlon,jlat
            CALL flush(6)
            WRITE(6,*)'JFRAC   PHGROUND   ZHGROUND   PTABGW_H(JFRAC)'
            WRITE(6,*)jfrac,phground(jlon,jlat),zhground(jlon,jlat),ptabgw_h(jlon,jlat,jfrac)
            CALL flush(6)
            CALL abort_trip('GWF_CPL_UPDATE:Problem with IINF or ISUP')
          ENDIF        
        ENDIF
     ENDDO
  ENDDO
!
  DO jlat=1,ilat
     DO jlon=1,ilon      
        IF(lmask(jlon,jlat))THEN
          pfwtd(jlon,jlat) = ptabgw_f(jlon,jlat,iinf(jlon,jlat))                                      &
                           + (zhground(jlon,jlat                )-ptabgw_h(jlon,jlat,iinf(jlon,jlat))) &
                           * (ptabgw_f(jlon,jlat,isup(jlon,jlat))-ptabgw_f(jlon,jlat,iinf(jlon,jlat))) &
                           / (ptabgw_h(jlon,jlat,isup(jlon,jlat))-ptabgw_h(jlon,jlat,iinf(jlon,jlat)))
          pfwtd(jlon,jlat) = min(1.0,pfwtd(jlon,jlat))
        ENDIF
     ENDDO
  ENDDO
!
ENDIF
!
!-------------------------------------------------------------------------------
! * Update the old groundwater height
!-------------------------------------------------------------------------------
!
phg_old(:,:)=phground(:,:)
!
IF (lhook) CALL dr_hook('GWF_CPL_UPDATE',1,zhook_handle)
!
END SUBROUTINE gwf_cpl_update