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(T, PTSTEP, 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 modd_teb_n, ONLY : teb_t
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_layer_e_budget_get_coef
!
IMPLICIT NONE
!
!
!*      0.1    Declarations of arguments
!
TYPE(teb_t), INTENT(INOUT) :: T
!
REAL              ,   INTENT(IN)  :: PTSTEP      ! time step
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
INTEGER :: ILAYER
!
REAL, DIMENSION(SIZE(PTDEEP_A),SIZE(T%XT_ROOF,2)) :: 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),SIZE(T%XT_ROOF,2)) :: ZW   ! work array
REAL, DIMENSION(SIZE(PTDEEP_A),SIZE(T%XT_ROOF,2)) :: 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.
!              ---------------------------------------------------------------
!
ilayer = SIZE(t%XT_ROOF,2)
!
 CALL layer_e_budget_get_coef(t%XT_ROOF, ptstep, zimpl, t%XHC_ROOF, t%XTC_ROOF, t%XD_ROOF, &
                              za, zb, zc, zy )
!
!-------------------------------------------------------------------------------
!
!*      2.0    Solving of the equation from bottom to top
!              ------------------------------------------
!
! layer at bottom of roof
!
zdet(:)   = zb(:,ilayer)
!
zt(:,ilayer) = zy(:,ilayer) / zdet(:)
!
! internal layers & top layer (but without the external heat flux term)
!
DO jk=ilayer-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 = T%XT_ROOF(:,1)
!PTDEEP_A = 0.
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('ROOF_IMPL_COEF',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE roof_impl_coef