SURFEX v7.3
General documentation of Surfex
|
00001 ! ############################################### 00002 SUBROUTINE ROOF_IMPL_COEF(PTSTEP,PTDEEP_A,PTDEEP_B) 00003 ! ############################################### 00004 ! 00005 !! 00006 !! PURPOSE 00007 !! ------- 00008 ! 00009 ! Computes the corfficients for implicitation of upper 00010 ! roof layer with what is above 00011 ! 00012 ! 00013 !!** METHOD 00014 ! ------ 00015 ! 00016 ! One computes a guess assuming a zero flux condition at the base 00017 ! of the roof. One solves the half part of the tridiagonal matrix 00018 ! fromm bottom to top. 00019 ! 00020 !! The classical tridiagonal algorithm is used to invert the 00021 !! implicit operator (from bottom to top only). Its matrix is given by: 00022 !! 00023 !! ( b(1) c(1) 0 0 0 0 0 0 ) 00024 !! ( a(2) b(2) c(2) 0 ... 0 0 0 0 ) 00025 !! ( 0 a(3) b(3) c(3) 0 0 0 0 ) 00026 !! ....................................................................... 00027 !! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 ) 00028 !! ....................................................................... 00029 !! ( 0 0 0 0 0 ... a(n-1) b(n-1) c(n-1)) 00030 !! ( 0 0 0 0 0 ... 0 a(n) b(n) ) 00031 !! 00032 !! 00033 !! AUTHOR 00034 !! ------ 00035 !! 00036 !! V. Masson * Meteo-France * 00037 !! 00038 !! MODIFICATIONS 00039 !! ------------- 00040 !! Original 01/2013 00041 !! 00042 !------------------------------------------------------------------------------- 00043 ! 00044 !* 0. DECLARATIONS 00045 ! ------------ 00046 ! 00047 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00048 USE PARKIND1 ,ONLY : JPRB 00049 ! 00050 USE MODD_TEB_n, ONLY : XD_ROOF, XTC_ROOF, XHC_ROOF, XT_ROOF, NROOF_LAYER 00051 ! 00052 USE MODI_LAYER_E_BUDGET_GET_COEF 00053 ! 00054 IMPLICIT NONE 00055 ! 00056 ! 00057 !* 0.1 Declarations of arguments 00058 ! 00059 REAL , INTENT(IN) :: PTSTEP ! time step 00060 REAL, DIMENSION(:), INTENT(OUT) :: PTDEEP_A, PTDEEP_B 00061 ! Deep soil temperature (prescribed) 00062 ! PTDEEP_A = Deep soil temperature 00063 ! coefficient depending on flux 00064 ! PTDEEP_B = Deep soil temperature (prescribed) 00065 ! which models heating/cooling from 00066 ! below the diurnal wave penetration 00067 ! (surface temperature) depth. If it 00068 ! is FLAGGED as undefined, then the zero 00069 ! flux lower BC is applied. 00070 ! Tdeep = PTDEEP_B + PTDEEP_A * PDEEP_FLUX 00071 ! (with PDEEP_FLUX in W/m2) 00072 ! 00073 !* 0.2 Local variables 00074 ! 00075 REAL :: ZIMPL = 1.0 ! implicit coefficient 00076 INTEGER :: JK ! loop counter 00077 ! 00078 REAL, DIMENSION(SIZE(PTDEEP_A),NROOF_LAYER) :: ZA, ! lower diag. 00079 ZB, ! main diag. 00080 ZC, ! upper diag. 00081 ZY ! r.h.s. 00082 00083 REAL, DIMENSION(SIZE(PTDEEP_A)) :: ZDET ! work array 00084 REAL, DIMENSION(SIZE(PTDEEP_A),NROOF_LAYER) :: ZW ! work array 00085 REAL, DIMENSION(SIZE(PTDEEP_A),NROOF_LAYER) :: ZT ! guess of T 00086 ! 00087 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00088 !------------------------------------------------------------------------------- 00089 IF (LHOOK) CALL DR_HOOK('ROOF_IMPL_COEF',0,ZHOOK_HANDLE) 00090 ! 00091 !* 1.0 Coefficients of the tridioagonal matrix for heat conduction eq. 00092 ! --------------------------------------------------------------- 00093 ! 00094 CALL LAYER_E_BUDGET_GET_COEF( XT_ROOF, PTSTEP, ZIMPL, XHC_ROOF, XTC_ROOF, XD_ROOF, & 00095 ZA, ZB, ZC, ZY ) 00096 ! 00097 !------------------------------------------------------------------------------- 00098 ! 00099 !* 2.0 Solving of the equation from bottom to top 00100 ! ------------------------------------------ 00101 ! 00102 ! layer at bottom of roof 00103 ! 00104 ZDET(:) = ZB(:,NROOF_LAYER) 00105 ! 00106 ZT (:,NROOF_LAYER) = ZY(:,NROOF_LAYER) / ZDET(:) 00107 ! 00108 ! internal layers & top layer (but without the external heat flux term) 00109 ! 00110 DO JK=NROOF_LAYER-1,1,-1 00111 ZW (:,JK) = ZA(:,JK+1)/ZDET(:) 00112 ZDET(:) = ZB(:,JK ) - ZC(:,JK)*ZW(:,JK) 00113 ZT (:,JK) = ( ZY(:,JK) - ZC(:,JK)*ZT(:,JK+1) ) / ZDET(:) ! + FLUX / ZDET 00114 ! ! for layer 1 00115 ! ! because the external 00116 ! ! flux would be 00117 ! ! included in the Y 00118 ! ! term 00119 END DO 00120 ! 00121 ! Implicit coefficients for the heat flux 00122 ! 00123 PTDEEP_B = ZT (:,1) 00124 PTDEEP_A = 1. / ZDET(:) 00125 ! 00126 !* The following lines are here if you want to test the explicit coupling 00127 !PTDEEP_B = XT_ROOF(:,1) 00128 !PTDEEP_A = 0. 00129 !------------------------------------------------------------------------------- 00130 IF (LHOOK) CALL DR_HOOK('ROOF_IMPL_COEF',1,ZHOOK_HANDLE) 00131 !------------------------------------------------------------------------------- 00132 ! 00133 END SUBROUTINE ROOF_IMPL_COEF