roof_impl_coef.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 roof_impl_coef(PTSTEP, KROOF_LAYER, PD_ROOF, PTC_ROOF, PHC_ROOF, PT_ROOF, PTDEEP_A,PTDEEP_B)
!   ###############################################
!
!!
!!    PURPOSE
!!    -------
!
!     Computes the corfficients for implicitation of upper
!     roof layer with what is above
!         
!     
!!**  METHOD
!     ------
!
!  One computes a guess assuming a zero flux condition at the base
!  of the roof. One solves the half part of the tridiagonal matrix
!  fromm bottom to top.
!
!!     The classical tridiagonal algorithm is used to invert the 
!!     implicit operator (from bottom to top only). Its matrix is given by:
!!
!!     (  b(1)      c(1)      0        0        0         0        0        0  )
!!     (  a(2)      b(2)     c(2)      0  ...    0        0        0        0  ) 
!!     (   0        a(3)     b(3)     c(3)       0        0        0        0  ) 
!!      .......................................................................
!!     (   0   ...   0      a(k)      b(k)     c(k)       0   ...  0        0  ) 
!!      .......................................................................
!!     (   0         0        0        0        0 ...  a(n-1)   b(n-1)   c(n-1))
!!     (   0         0        0        0        0 ...     0      a(n)     b(n) )
!!
!!
!!    AUTHOR
!!    ------
!!
!!      V. Masson           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2013
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_layer_e_budget_get_coef
!
IMPLICIT NONE
!
!
!*      0.1    Declarations of arguments
!
REAL              ,   INTENT(IN)  :: ptstep      ! time step
INTEGER           ,   INTENT(IN)  :: kroof_layer ! number of roof layers
REAL, DIMENSION(:,:), INTENT(IN)  :: pd_roof   ! thickness            of each layer
REAL, DIMENSION(:,:), INTENT(IN)  :: ptc_roof  ! thermal conductivity of each layer
REAL, DIMENSION(:,:), INTENT(IN)  :: phc_roof  ! heat capacity        of each layer
REAL, DIMENSION(:,:), INTENT(IN)  :: pt_roof   ! temperature          of each layer
REAL, DIMENSION(:),   INTENT(OUT) :: ptdeep_a, ptdeep_b     
                                              ! Deep soil temperature (prescribed)
!                                      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)
!
!*      0.2    Local variables
!
REAL :: zimpl = 1.0        ! implicit coefficient
INTEGER :: jk              ! loop counter
!
REAL, DIMENSION(SIZE(PTDEEP_A),KROOF_LAYER) :: za,& ! lower diag.
                                               zb,& ! main  diag.
                                               zc,& ! upper diag.
                                               zy   ! r.h.s.

REAL, DIMENSION(SIZE(PTDEEP_A))             :: zdet ! work array
REAL, DIMENSION(SIZE(PTDEEP_A),KROOF_LAYER) :: zw   ! work array
REAL, DIMENSION(SIZE(PTDEEP_A),KROOF_LAYER) :: zt   ! guess of T
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('ROOF_IMPL_COEF',0,zhook_handle)
!
!*      1.0    Coefficients of the tridioagonal matrix for heat conduction eq.
!              ---------------------------------------------------------------
!
 CALL layer_e_budget_get_coef( pt_roof, ptstep, zimpl, phc_roof, ptc_roof, pd_roof, &
                              za, zb, zc, zy )
!
!-------------------------------------------------------------------------------
!
!*      2.0    Solving of the equation from bottom to top
!              ------------------------------------------
!
! layer at bottom of roof
!
zdet(:)   = zb(:,kroof_layer)
!
zt(:,kroof_layer) = zy(:,kroof_layer) / zdet(:)
!
! internal layers & top layer (but without the external heat flux term)
!
DO jk=kroof_layer-1,1,-1
  zw(:,jk)  = za(:,jk+1)/zdet(:)
  zdet(:)     = zb(:,jk  ) - zc(:,jk)*zw(:,jk)
  zt(:,jk)  = ( zy(:,jk) - zc(:,jk)*zt(:,jk+1) ) / zdet(:) ! + FLUX / ZDET
  !                                                          ! for layer 1
  !                                                          ! because the external
  !                                                          ! flux would be
  !                                                          ! included in the Y
  !                                                          ! term
END DO
!
! Implicit coefficients for the heat flux
!
ptdeep_b = zt(:,1)
ptdeep_a = 1. / zdet(:) 
!
!* The following lines are here if you want to test the explicit coupling
!PTDEEP_B = PT_ROOF(:,1)
!PTDEEP_A = 0.
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('ROOF_IMPL_COEF',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE roof_impl_coef