mode_hydro_dif.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.
!     ######spl
      MODULE mode_hydro_dif 
!     ################
!
!!****  *MODE_HYDRO_DIF * - pedo-transfert functions
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!    
!!
!!    EXTERNAL
!!    --------
!!       
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!      B. Decharme       * Meteo France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original        11/2010
!-----------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!
!
INTERFACE vapcondcf
  MODULE PROCEDURE vapcondcf       
END INTERFACE
!
INTERFACE infmax_func
  MODULE PROCEDURE infmax_func
END INTERFACE
!
INTERFACE tridiag_dif
  MODULE PROCEDURE tridiag_dif
END INTERFACE
!
!-------------------------------------------------------------------------------
CONTAINS
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!vapor conductivity (m s-1)
!-------------------------------------------------------------------------------
!
FUNCTION vapcondcf(PTG,PPS,PWG,PWGI,PPSIA,PWSAT,PWFC,PQSAT,PQSATI,KWG_LAYER,KNL) RESULT(PVAPCOND)
!
! Uses method of Braud et al. (1993) for
!
USE modd_csts,       ONLY : xmv, xmd, xtt, xp00, xg, xrv, xrholw
USE modd_isba_par,   ONLY : xwgmin
USE modd_surf_par,   ONLY : xundef
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:  ), INTENT(IN)         :: PPS
REAL, DIMENSION(:,:), INTENT(IN)         :: PWG,PWGI,PPSIA,PWSAT, &
                                            PWFC,PTG,PQSAT,PQSATI
INTEGER, DIMENSION(:), INTENT(IN)        :: KWG_LAYER       !Moisture layer
INTEGER,               INTENT(IN)        :: KNL             ! number of vertical levels
!
REAL, DIMENSION(SIZE(PWG,1),SIZE(PWG,2)) :: PVAPCOND
!
!
!*      0.2    declarations of local variables
!
REAL    :: ZDVA, ZFVA, ZCHI, ZHUM, ZWORK,  &
           ZPV, ZESAT, ZESATI, ZWG, ZVC
!
INTEGER :: INI, JJ, JL, IDEPTH
!
! Parameters:
!
REAL, PARAMETER                     :: ZTORTY = 0.66         ! (-)
REAL, PARAMETER                     :: ZNV    = 1.88         ! (-)
REAL, PARAMETER                     :: ZCV    = 2.17e-5      ! (m2/s)
REAL, PARAMETER                     :: ZWK    = 0.05         ! (m3 m-3)
REAL, PARAMETER                     :: ZLIM   = tiny(1.0)
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_HYDRO_DIF:VAPCONDCF',0,zhook_handle)
!
ini = SIZE(pwg,1)
!
pvapcond(:,:) = 0.0
!
! Only perform this computation if the soil is sufficiently
! dry (as otherwise the hydraulic conductivity dominates
! the diffusion coefficient). Arbitrarily base threshold on field
! capacity water content:
!
DO jl=1,knl
   DO jj=1,ini
!
      idepth = kwg_layer(jj)
      zwg    = pwg(jj,jl) + pwgi(jj,jl)
!      
      IF(jl<=idepth .AND. zwg < pwfc(jj,jl) .AND. zwg > xwgmin)THEN
!
!        Vapor pressure over liquid and solid water surfaces (Pa), respectively:
!
         zesat  = pqsat(jj,jl)* pps(jj)/((xmv/xmd)+pqsat(jj,jl) *(1.-(xmv/xmd)))
!
         zesati = pqsati(jj,jl)*pps(jj)/((xmv/xmd)+pqsati(jj,jl)*(1.-(xmv/xmd)))
!
!        molecular diffusivity of water vapor (m2 s-1):
!
         zwork  = znv*log(ptg(jj,jl)/xtt)
         zdva   = zcv*(xp00/pps(jj))*exp(zwork)
!
!        function of pore space: 
!
         zfva   = (pwsat(jj,jl) - zwg)*(1.+(zwg/(pwsat(jj,jl)-zwk)))
         zfva   = min(zfva,pwsat(jj,jl))
!
!        relative humidity of air in soil pores:
!
         zhum   = max(zlim,exp(ppsia(jj,jl)*xg/(xrv*ptg(jj,jl))))
!
!        fraction of frozen water:
!
         zchi   = pwgi(jj,jl)/zwg
!
!        vapor pressure within pore space (Pa):
!
         zpv    = zhum*(zchi*zesat + (1.-zchi)*zesati)
!
!        vapor conductivity (kg m-2 s-1)
!
         zvc    = ztorty*pps(jj)*zdva*zfva*xg*zpv/                  &
                  ((pps(jj)-zpv)*(xrv*xrv*ptg(jj,jl)*ptg(jj,jl)))
!
!        vapor conductivity (m s-1)
!
         pvapcond(jj,jl) = zvc/xrholw
!
      ENDIF
!
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('MODE_HYDRO_DIF:VAPCONDCF',1,zhook_handle)
!
END FUNCTION vapcondcf
!
!-------------------------------------------------------------------------------
!Green-Ampt approximation (derived form) for maximum infiltration
!-------------------------------------------------------------------------------
!
FUNCTION infmax_func(PWG,PWSAT,PFRZ,PCONDSAT,PMPOTSAT,PBCOEF,PDZG,PDG,KLAYER_HORT)
USE yomhook      ,ONLY : lhook,   dr_hook
USE modd_sgh_par, ONLY : xhort_depth
USE parkind1     ,ONLY : jprb
IMPLICIT NONE
REAL, DIMENSION(:,:), INTENT(IN) :: PWG
REAL, DIMENSION(:,:), INTENT(IN) :: PWSAT           
REAL, DIMENSION(:,:), INTENT(IN) :: PFRZ
REAL, DIMENSION(:,:), INTENT(IN) :: PCONDSAT            
REAL, DIMENSION(:,:), INTENT(IN) :: PMPOTSAT    
REAL, DIMENSION(:,:), INTENT(IN) :: PBCOEF 
REAL, DIMENSION(:,:), INTENT(IN) :: PDZG
REAL, DIMENSION(:,:), INTENT(IN) :: PDG
INTEGER,              INTENT(IN) :: KLAYER_HORT   
!
REAL, DIMENSION(SIZE(PWG,1)) :: ZGREEN_AMPT, ZDEPTH
REAL                         :: ZS, ZCOEF
INTEGER                      :: JJ,JL,INI
!
REAL, DIMENSION(SIZE(PWG,1)) :: INFMAX_FUNC
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('MODE_HYDRO_DIF:INFMAX_FUNC',0,zhook_handle)
!
ini   =SIZE(pwg,1)
!
zgreen_ampt(:) = 0.0
zdepth(:) = 0.0
!
DO jl=1,klayer_hort
   DO jj=1,ini  
      IF(zdepth(jj)<xhort_depth)THEN
         zs              = min(1.0,pwg(jj,jl)/pwsat(jj,jl))
         zcoef           = pbcoef(jj,jl)*pmpotsat(jj,jl)*(zs-1.0)/pdzg(jj,jl)        
         zgreen_ampt(jj) = zgreen_ampt(jj)+pdzg(jj,jl)*pfrz(jj,jl)*pcondsat(jj,jl)*(zcoef+1.0)
         zdepth(jj) = pdg(jj,jl)
      ENDIF      
   ENDDO
ENDDO
!
infmax_func(:) = zgreen_ampt(:)/zdepth(:)
!
IF (lhook) CALL dr_hook('MODE_HYDRO_DIF:INFMAX_FUNC',1,zhook_handle)
END FUNCTION infmax_func
!
!-------------------------------------------------------------------------------
!Solve tridiagonal matrix (for method see tridiag_ground.F90)
!-------------------------------------------------------------------------------
!
SUBROUTINE tridiag_dif(PAMTRX,PBMTRX,PCMTRX,PFRC,KWG_LAYER,KNL,PSOL)
USE modd_surf_par, ONLY : xundef
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
IMPLICIT NONE
REAL,    DIMENSION(:,:), INTENT(IN)  :: PAMTRX    ! lower diag. elements of A matrix
REAL,    DIMENSION(:,:), INTENT(IN)  :: PBMTRX    ! main  diag. elements of A matrix
REAL,    DIMENSION(:,:), INTENT(IN)  :: PCMTRX    ! upper diag. elements of A matrix
REAL,    DIMENSION(:,:), INTENT(IN)  :: PFRC      ! Forcing term
INTEGER, DIMENSION(:),   INTENT(IN)  :: KWG_LAYER !Moisture layer
INTEGER,                 INTENT(IN)  :: KNL       ! number of vertical levels
REAL,    DIMENSION(:,:), INTENT(OUT) :: PSOL      ! solution of A.SOL = FRC
!
REAL, DIMENSION(SIZE(PFRC,1),SIZE(PFRC,2)) :: ZWORK! work array
REAL, DIMENSION(SIZE(PFRC,1))              :: ZDET ! work array
INTEGER                                    :: JL   ! vertical loop control
INTEGER :: JJ, INI, IDEPTH
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('MODE_HYDRO_DIF:TRIDIAG_DIF',0,zhook_handle)
!
ini=SIZE(pfrc,1)
!
psol(:,:)=xundef
zwork(:,:)=xundef
!
!first level
zdet(:)   = pbmtrx(:,1)
psol(:,1) = pfrc(:,1) / zdet(:)
!
!other levels
DO jl=2,knl
   DO jj=1,ini
      idepth=kwg_layer(jj)
      IF(jl<=idepth)THEN
        zwork(jj,jl) = pcmtrx(jj,jl-1)/zdet(jj)
        zdet(jj)    = pbmtrx(jj,jl) - pamtrx(jj,jl)*zwork(jj,jl)
        psol(jj,jl) = (pfrc(jj,jl) - pamtrx(jj,jl)*psol(jj,jl-1))/zdet(jj)  
      ENDIF
   ENDDO 
ENDDO        
!
!levels going down
DO jl=knl-1,1,-1
   DO jj=1,ini
      idepth=kwg_layer(jj)
      IF(jl<idepth)THEN
         psol(jj,jl) = psol(jj,jl)-zwork(jj,jl+1)*psol(jj,jl+1)
      ENDIF
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('MODE_HYDRO_DIF:TRIDIAG_DIF',1,zhook_handle)
!
END SUBROUTINE tridiag_dif
!
!-------------------------------------------------------------------------------
!
END MODULE mode_hydro_dif