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