tridiag_ground_rm_coefs.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.
!     #########
SUBROUTINE tridiag_ground_rm_coefs(PTSTEP,PDEPTH,PTEMP,PHEATCAP,PCONDTRM,              &
     psource,ptdeep_a,ptdeep_b,pconda_delz,pa_coef,pb_coef)
!
!
!!****  *TRIDIAG_GROUND_RM_COEF*  
!!
!!    PURPOSE
!!    -------
!
!     This routine computes the coefficients of a Tridiagnoal matrix using
!     the method of Richtmeyer and Morton (1967): this routine corresponds to the **upward sweep**
!     The lower boundary condition is based on the linear eq ==> TN = B + A FN
!     Where B and A are imposed and correspond to either a Dirichlet (prescribed lower boundary T)
!     or a lower boundary flux (Neumann type). If B is FLAGGED, then
!     lower boundary flux is imposed as zero.
!     This routine is used to
!     eliminate T2(t+dt) (sub-surface T at time t+dt) from an energy balance Eq (where T1=Tsfc)
!     NOTE the solution is computed from TRIDIAG_GROUND_RM_SOLN.F90
!         
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!
!     Richtmeyer, R. and K. Morton, 1967: Difference method for initial values problems,
!     Interscience Publishers, 2.
!
!!    AUTHOR
!!    ------
!!
!!      A. Boone           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    21/03/11 
!!      Modif       23/02/12 A. Boone: Optimization
!!      Modif       03/2013  A. Boone: MEB
!!      Modif       06/2013  A. Boone: use TN=B+A FN form for lower BC
!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_surf_par,   ONLY : xundef
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    Declarations of arguments
!
REAL,                 INTENT(IN)  :: ptstep        ! time-step                         (s)
REAL, DIMENSION(:,:), INTENT(IN)  :: pdepth        ! soil layer depth                  (m)
REAL, DIMENSION(:,:), INTENT(IN)  :: ptemp         ! surface and sub-surface soil 
!                                                  ! temperature profile               (K)
REAL, DIMENSION(:,:), INTENT(IN)  :: pheatcap      ! soil heat capacity                (J/K/m3)
REAL, DIMENSION(:,:), INTENT(IN)  :: pcondtrm      ! soil thermal conductivity         (W/m/K)
REAL, DIMENSION(:,:), INTENT(IN)  :: psource       ! soil heat source function         (J/m2/s)
REAL, DIMENSION(:),   INTENT(IN)  :: ptdeep_a, ptdeep_b
!                                                  ! PTDEEP_A = Deep soil temperature
!                                                  ! coefficient depending on flux     (m2/K/W)
!                                                  ! PTDEEP_B = Deep soil temperature  (K)
!                                                  ! (prescribed)
!                                                  ! which models heating/cooling from
!                                                  ! below the diurnal wave penetration
!                                                  ! (surface temperature) depth. If it
!                                                  ! is FLAGGED as undefined, then the zero
!                                                  ! flux lower BC is applied.
!                                                  ! Tdeep = PTDEEP_B + PTDEEP_A * PDEEP_FLUX
!                                                  ! (with PDEEP_FLUX in W/m2)
!
REAL, DIMENSION(:),   INTENT(OUT) :: pconda_delz   ! ratio: ground flux thermal 
                                                   ! conductivity/ 
                                                   ! ground flux thickness             (W/m2/K)
REAL, DIMENSION(:,:), INTENT(OUT) :: pa_coef       ! RM67 A-soil coefficient           (-)
REAL, DIMENSION(:,:), INTENT(OUT) :: pb_coef       ! RM67 B-soil coefficient           (K)
!
!
!*      0.2    Declarations of local variables
!
INTEGER                                        :: jj, ji, inl, ini
REAL, DIMENSION(SIZE(PDEPTH,1),SIZE(PDEPTH,2)) :: zd_g             ! soil layer thicknesses             (m)
REAL, DIMENSION(SIZE(PDEPTH,1),SIZE(PDEPTH,2)) :: zdlz             ! thickness between layer mid_points (m)
REAL, DIMENSION(SIZE(PDEPTH,1),SIZE(PDEPTH,2)) :: zlambda          ! ratio of thermal cond to dz        (W/m2/K)
REAL, DIMENSION(SIZE(PDEPTH,1),SIZE(PDEPTH,2)) :: zgheatcap        ! effective heat capacity         
!                                                                  ! coefficient                        (J/K/m2)
REAL, DIMENSION(SIZE(PDEPTH,1),SIZE(PDEPTH,2)) :: zthrm            ! interfacial therm. cond.           (W/m/K) 
REAL, DIMENSION(SIZE(PDEPTH,1))                :: zc_coef          ! working denominator for coefs
REAL, DIMENSION(SIZE(PDEPTH,1))                :: za_coefd         !
REAL, DIMENSION(SIZE(PDEPTH,1))                :: zb_coefd         ! 
REAL, DIMENSION(SIZE(PDEPTH,1),SIZE(PDEPTH,2)) :: zsink            ! sink term (can be input)
!
REAL(KIND=JPRB) :: zhook_handle
!
!*      0.3    Declarations of local parameters
!
REAL, PARAMETER                                :: zdz_min = 1.e-6  ! m
!
!---------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TRIDIAG_GROUND_RM_COEFS',0,zhook_handle)
!
!*       0.     Initialization
!               --------------
!
ini              = SIZE(pdepth,1)
inl              = SIZE(pdepth,2)
!
!*       1.     Compute needed parameters
!               -------------------------
! Needed parameters: interfacial thermal conductivity,
! layer thicknesses, and layer average heat capacities:

! Soil layer thicknesses, dz, from depths (base of each layer):

zd_g(:,1)         = pdepth(:,1)
DO jj=2,inl
   DO ji=1,ini
      zd_g(ji,jj) = pdepth(ji,jj) - pdepth(ji,jj-1)  
   ENDDO
ENDDO
zd_g(:,:)         = max(zdz_min, zd_g(:,:)) ! just for numerical reasons

! distance between layer mid_points (m):

DO jj=1,inl-1
   DO ji=1,ini
      zdlz(ji,jj)  = (zd_g(ji,jj)+zd_g(ji,jj+1))*0.5
   ENDDO
ENDDO
zdlz(:,inl)       = zd_g(:,inl)*0.5

! effective heat capacity coefficient (J/K/m2):

zgheatcap(:,:)    = zd_g(:,:)*pheatcap(:,:)

! interfacial thermal conductivity (W/m/K):

DO jj=1,inl-1
   DO ji=1,ini
      zthrm(ji,jj)  = (zd_g(ji,jj)+zd_g(ji,jj+1))/              &
                     ((zd_g(ji,jj+1)/pcondtrm(ji,jj+1))+        &
                     (zd_g(ji,jj)  /pcondtrm(ji,jj)  ))
   ENDDO
ENDDO
zthrm(:,inl)        = pcondtrm(:,inl) 
!
!
! Energy sink: 
! NOTE for now, set this part to zero as accounted for elsewhere.
!!!ZSINK(:,:)       = -PPHASE(:,:)*PTSTEP/ZGHEATCAP(:,:)
!
zsink(:,:)          = -psource(:,:)  ! J/m2/s
!
! Thermal cond/dz ratio (W m-2 K-1):
!
zlambda(:,:)        = zthrm(:,:)/zdlz(:,:)
!
! Save sfc interfacial thermal cond/dz ratio (W m-2 K-1):
!
pconda_delz(:)      = zlambda(:,1)
!
!
!*       2.     COEFFICIENTS
!               ------------
!
! Start at base of soil or snow etc...assuming an imposed lower boundary flux,
! and move up towards surface:

! initialize:

pa_coef(:,:)        = 0.0
pb_coef(:,:)        = 0.0

! lowest layer: flux is prescribed (can be 0)
! where GBFLUX = tcond/dz ( T_N - T_N+1 )  (W m-2)

WHERE(ptdeep_b(:) == xundef) ! no flux at model base

   zc_coef(:)     = (zgheatcap(:,inl)/ptstep)                           &
                   + zlambda(:,inl-1)
   pa_coef(:,inl) = zlambda(:,inl-1)/zc_coef(:)
   pb_coef(:,inl) = ( ptemp(:,inl)*(zgheatcap(:,inl)/ptstep)            &
                   - zsink(:,inl) )/zc_coef(:)

   za_coefd(:)    = xundef
   zb_coefd(:)    = xundef
   
ELSEWHERE 

! corresponds to Dirichlet or Neumann type lower BC (non-zero flux)

   zc_coef(:)     = 1. - ptdeep_a(:)*zlambda(:,inl)
   za_coefd(:)    = -zlambda(:,inl)/zc_coef(:)
   zb_coefd(:)    =  zlambda(:,inl)*ptdeep_b(:)/zc_coef(:)

   zc_coef(:)     = (zgheatcap(:,inl)/ptstep)                          &
                    + zlambda(:,inl-1) - za_coefd(:)
   pa_coef(:,inl) = zlambda(:,inl-1)/zc_coef(:)
   pb_coef(:,inl) = ( ptemp(:,inl)*(zgheatcap(:,inl)/ptstep)           &
                    + zb_coefd(:) - zsink(:,inl) )/zc_coef(:)

END WHERE

! interior layers:

DO jj=inl-1,2,-1
   DO ji=1,ini
      zc_coef(ji)    = (zgheatcap(ji,jj)/ptstep)                       &
                      + zlambda(ji,jj)*(1.0-pa_coef(ji,jj+1))          &
                      + zlambda(ji,jj-1)
      pa_coef(ji,jj) = zlambda(ji,jj-1)/zc_coef(ji)
      pb_coef(ji,jj) = ( ptemp(ji,jj)*(zgheatcap(ji,jj)/ptstep)        &
                      + zlambda(ji,jj)*pb_coef(ji,jj+1)                &
                      - zsink(ji,jj) )/zc_coef(ji)
   ENDDO
ENDDO
!
! NOTE: uppermost coefficients have been kept at
!       zero (initial value)/not computed as they're 
!       not used in computations (implicit in surface E budget
!       computed elsewhere)
!
IF (lhook) CALL dr_hook('TRIDIAG_GROUND_RM_COEFS',1,zhook_handle)
!-------------------------------------------------------------------
END SUBROUTINE tridiag_ground_rm_coefs