soil_heatdif.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 soil_heatdif(PTSTEP,PDZG,PDZDIF,PSOILCONDZ,   &
                              PSOILHCAPZ,PCT,PTERM1,PTERM2,    &
                              PTDEEP_A,PTDEEP_B,PTG,PDEEP_FLUX,&
                              PFLUX_COR                        )  
!     ############################################################################
!
!!****  *SOIL_HEATDIF*  
!
!!    PURPOSE
!!    -------
!     This subroutine solves the 1-d diffusion of 'PTG' using a
!     layer averaged set of equations which are time differenced
!     using the backward-difference scheme (implicit). 
!     The eqs are solved rapidly by taking advantage of the
!     fact that the matrix is tridiagonal. This is a very
!     general routine and can be used for the 1-d diffusion of any 
!     quantity as long as the diffusity is not a function of the
!     quantity being diffused. Aaron Boone 8-98. Soln to the eq:
!
!                    dQ    d    dQ       
!                 c  -- =  -- k -- 
!                    dt    dx   dx    
!
!     where k = k(x) (thermal conductivity), c = c(x) (heat capacity)
!     Diffusivity is k/c
!
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    USE MODD_PARAMETERS
!!    USE MODI_TRIDIAG_GROUND
!!      
!!    REFERENCE
!!    ---------
!!
!!    Boone (2000)
!!    Boone et al. (2000)
!!      
!!    AUTHOR
!!    ------
!!      A. Boone          * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    16/02/00   Boone
!!      Modif       08/2011 B. Decharme : Optimization
!!      Modif       10/2014 B. Decharme : Use harmonic mean to compute 
!!                                        interfacial thermal conductivities
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_surf_par,   ONLY : xundef
!
USE modi_tridiag_ground
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!
REAL, INTENT(IN)                    :: PTSTEP ! Model time step (s)
!
REAL, DIMENSION(:), INTENT(IN)      :: PCT, PTERM1, PTERM2, PTDEEP_A, PTDEEP_B
!                                      PCT    = thermal inertia [(m2 K)/J]
!                                      PTERM1 = coefficient of linearization
!                                               of surface energy budget 
!                                      PTERM2 = coefficient of linearization
!                                               of surface energy budget 
!                                      PTDEEP_A = Deep soil temperature
!                                                 coefficient depending on flux
!                                      PTDEEP_B = Deep soil temperature (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(IN)    :: PSOILCONDZ, PSOILHCAPZ, PFLUX_COR
!                                      PSOILCONDZ = soil thermal conductivity [W/(K m)]
!                                      PSOILHCAPZ = soil heat capacity [J/(m3 K)]
!                                      PFLUX_COR  = correction flux to close energy budget (W/m2)
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PDZDIF, PDZG
!                                      PDZDIF = distance between consecuative layer mid-points
!                                      PDZG   = soil layers thicknesses
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PTG
!                                      PTG    = soil temperature (K)
!
REAL, DIMENSION(:), INTENT(OUT)     :: PDEEP_FLUX ! Heat flux at bottom of ISBA (W/m2)
!
!
!*      0.2    declarations of local variables
!
INTEGER                                  :: JJ, JL
!
INTEGER                                  :: INI, INLVLD ! Number of point and grid layers
!
REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2)) :: ZTGM, ZDTERM, ZCTERM,   &
                                                    ZFRCV, ZAMTRX, ZBMTRX,     &
                                                    ZCMTRX  
!
REAL :: ZWORK1, ZWORK2
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
! Initialize local variables:
!
IF (lhook) CALL dr_hook('SOIL_HEATDIF',0,zhook_handle)
zdterm(:,:)    = 0.0
zcterm(:,:)    = 0.0
zfrcv(:,:)     = 0.0
zamtrx(:,:)    = 0.0
zbmtrx(:,:)    = 0.0
zcmtrx(:,:)    = 0.0
ztgm(:,:)   = ptg(:,:)
!
ini    = SIZE(ptg(:,:),1)
inlvld = SIZE(ptg(:,:),2)
!
!-------------------------------------------------------------------------------
!
! Calculate tri-diagonal matrix coefficients:
!
DO jl=1,inlvld
  DO jj=1,ini
     IF(jl<inlvld)THEN
       zwork1 = pdzg(jj,jl  )/(2.0*pdzdif(jj,jl)*psoilcondz(jj,jl  ))
       zwork2 = pdzg(jj,jl+1)/(2.0*pdzdif(jj,jl)*psoilcondz(jj,jl+1))
     ELSE
       zwork1 = pdzg(jj,jl)/(2.0*pdzdif(jj,jl)*psoilcondz(jj,jl))
       zwork2 = 0.0
     ENDIF
     zdterm(jj,jl)=1.0/(pdzdif(jj,jl)*(zwork1+zwork2))
     zcterm(jj,jl)= psoilhcapz(jj,jl)*pdzg(jj,jl)/ptstep
  ENDDO
ENDDO 
!
! - - -------------------------------------------------
!
! Upper BC
!
zamtrx(:,1) =  0.0
zbmtrx(:,1) =  1./(pct(:)*ptstep)
zcmtrx(:,1) = -pterm2(:)*zbmtrx(:,1)
zfrcv(:,1) =  pterm1(:)*zbmtrx(:,1)
!
! Interior Grid
!
DO jl=2,inlvld-1
   DO jj=1,ini
   zamtrx(jj,jl) = -zdterm(jj,jl-1) 
   zbmtrx(jj,jl) =  zcterm(jj,jl) + zdterm(jj,jl-1) + zdterm(jj,jl)
   zcmtrx(jj,jl) = -zdterm(jj,jl)
   zfrcv(jj,jl) =  zcterm(jj,jl)*ztgm(jj,jl)+pflux_cor(jj,jl)
   ENDDO
ENDDO
!
! Lower BC: 2 currently accounted for, Either zero flux
! or a fixed temperature 'TDEEP' 
!
zamtrx(:,inlvld) = -zdterm(:,inlvld-1) 
zcmtrx(:,inlvld) =  0.0                
!

!ZDEEP_FLUX=ZDTERM(:,INLVLD)*(ZTGM(:,INLVLD)-ZTDEEP(:))
!ZTDEEP=PTDEEP_A*ZDEEP_FLUX+PTDEEP_B
!
!ZDEEP_FLUX=ZDTERM(:,INLVLD)*(ZTGM(:,INLVLD)-PTDEEP_A*ZDEEP_FLUX-PTDEEP_B)
!ZDEEP_FLUX=ZDTERM(:,INLVLD)*(ZTGM(:,INLVLD)-PTDEEP_B)/(1.+ZDTERM(:,INLVLD)*PTDEEP_A)
!
!
WHERE(ptdeep_b(:) /= xundef)
!   ZBMTRX(:,INLVLD) =  ZCTERM(:,INLVLD) + ZDTERM(:,INLVLD-1) + ZDTERM(:,INLVLD)
!   ZFRCV(:,INLVLD)  =  ZCTERM(:,INLVLD)*ZTGM(:,INLVLD) + ZDTERM(:,INLVLD)*PTDEEP(:)
   zbmtrx(:,inlvld) =  zcterm(:,inlvld) + zdterm(:,inlvld-1) &
                     + zdterm(:,inlvld)/(1.+zdterm(:,inlvld)*ptdeep_a)
   zfrcv(:,inlvld)  =  zcterm(:,inlvld)*ztgm(:,inlvld) &
                     + zdterm(:,inlvld)*ptdeep_b(:)/(1.+zdterm(:,inlvld)*ptdeep_a) &
                     + pflux_cor(:,inlvld)
ELSEWHERE
   zbmtrx(:,inlvld) =  zcterm(:,inlvld) + zdterm(:,inlvld-1) 
   zfrcv(:,inlvld)  =  zcterm(:,inlvld)*ztgm(:,inlvld)+pflux_cor(:,inlvld)
END WHERE
!
! - - -------------------------------------------------
!
! Compute ZTGM (solution vector) 
! used for systems of equations involving tridiagonal 
! matricies.
!
 CALL tridiag_ground(zamtrx,zbmtrx,zcmtrx,zfrcv,ztgm)
!
!
! Update values in time:
!
ptg(:,:) = ztgm(:,:)
!
!* Deep soil Flux
!
pdeep_flux(:) = 0.
WHERE(ptdeep_b(:) /= xundef)
  pdeep_flux=zdterm(:,inlvld)*(ztgm(:,inlvld)-ptdeep_b)/(1.+zdterm(:,inlvld)*ptdeep_a)
END WHERE
!

IF (lhook) CALL dr_hook('SOIL_HEATDIF',1,zhook_handle)
!
!
!
END SUBROUTINE soil_heatdif