|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE TRIDIAG_GROUND(PA,PB,PC,PY,PX) 00003 ! ######################################### 00004 ! 00005 ! 00006 !!**** *TRIDIAG_GROUND* - routine to solve a time implicit scheme 00007 !! 00008 !! 00009 !! PURPOSE 00010 !! ------- 00011 ! The purpose of this routine is to resolve the linear system: 00012 ! 00013 ! A.X = Y 00014 ! 00015 ! where A is a tridiagonal matrix, and X and Y two vertical vectors. 00016 ! However, the computations are performed at the same time for all 00017 ! the verticals where an inversion of the system is necessary. 00018 ! This explain the dimansion of the input variables. 00019 ! 00020 !!** METHOD 00021 !! ------ 00022 !! 00023 !! Then, the classical tridiagonal algorithm is used to invert the 00024 !! implicit operator. Its matrix is given by: 00025 !! 00026 !! ( b(1) c(1) 0 0 0 0 0 0 ) 00027 !! ( a(2) b(2) c(2) 0 ... 0 0 0 0 ) 00028 !! ( 0 a(3) b(3) c(3) 0 0 0 0 ) 00029 !! ....................................................................... 00030 !! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 ) 00031 !! ....................................................................... 00032 !! ( 0 0 0 0 0 ... a(n-1) b(n-1) c(n-1)) 00033 !! ( 0 0 0 0 0 ... 0 a(n) b(n) ) 00034 !! 00035 !! 00036 !! All these computations are purely vertical and vectorizations are 00037 !! easely achieved by processing all the verticals in parallel. 00038 !! 00039 !! EXTERNAL 00040 !! -------- 00041 !! 00042 !! NONE 00043 !! 00044 !! IMPLICIT ARGUMENTS 00045 !! ------------------ 00046 !! 00047 !! REFERENCE 00048 !! --------- 00049 !! 00050 !! AUTHOR 00051 !! ------ 00052 !! V. Masson 00053 !! 00054 !! MODIFICATIONS 00055 !! ------------- 00056 !! Original May 13, 1998 00057 !! Modified : 00058 !! B. Decharme 08/12 Loop optimization 00059 !! --------------------------------------------------------------------- 00060 ! 00061 !* 0. DECLARATIONS 00062 ! 00063 ! 00064 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00065 USE PARKIND1 ,ONLY : JPRB 00066 ! 00067 IMPLICIT NONE 00068 ! 00069 ! 00070 !* 0.1 declarations of arguments 00071 ! 00072 REAL, DIMENSION(:,:), INTENT(IN) :: PA ! lower diag. elements of A matrix 00073 REAL, DIMENSION(:,:), INTENT(IN) :: PB ! main diag. elements of A matrix 00074 REAL, DIMENSION(:,:), INTENT(IN) :: PC ! upper diag. elements of A matrix 00075 REAL, DIMENSION(:,:), INTENT(IN) :: PY ! r.h.s. term 00076 ! 00077 REAL, DIMENSION(:,:), INTENT(OUT) :: PX ! solution of A.X = Y 00078 ! 00079 !* 0.2 declarations of local variables 00080 ! 00081 INTEGER :: JI ! number of point loop control 00082 INTEGER :: JK ! vertical loop control 00083 INTEGER :: INI ! number of point 00084 INTEGER :: INL ! number of vertical levels 00085 ! 00086 REAL, DIMENSION(SIZE(PA,1) ) :: ZDET ! work array 00087 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2)) :: ZW ! work array 00088 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00089 ! --------------------------------------------------------------------------- 00090 ! 00091 IF (LHOOK) CALL DR_HOOK('TRIDIAG_GROUND',0,ZHOOK_HANDLE) 00092 INI=SIZE(PX,1) 00093 INL=SIZE(PX,2) 00094 ! 00095 !* 1. levels going up 00096 ! --------------- 00097 ! 00098 !* 1.1 first level 00099 ! ----------- 00100 ! 00101 ZDET(:) = PB(:,1) 00102 00103 PX (:,1) = PY(:,1) / ZDET(:) 00104 ! 00105 !* 1.2 other levels 00106 ! ------------ 00107 ! 00108 DO JK=2,INL 00109 DO JI=1,INI 00110 ZW (JI,JK) = PC(JI,JK-1)/ZDET(JI) 00111 ZDET(JI) = PB(JI,JK ) - PA(JI,JK)*ZW(JI,JK) 00112 PX (JI,JK) = ( PY(JI,JK) - PA(JI,JK)*PX(JI,JK-1) ) / ZDET(JI) 00113 END DO 00114 END DO 00115 ! 00116 !------------------------------------------------------------------------------- 00117 ! 00118 !* 2. levels going down 00119 ! ----------------- 00120 ! 00121 DO JK=INL-1,1,-1 00122 DO JI=1,INI 00123 PX (JI,JK) = PX(JI,JK) - ZW(JI,JK+1)*PX(JI,JK+1) 00124 END DO 00125 END DO 00126 IF (LHOOK) CALL DR_HOOK('TRIDIAG_GROUND',1,ZHOOK_HANDLE) 00127 ! 00128 !------------------------------------------------------------------------------- 00129 ! 00130 END SUBROUTINE TRIDIAG_GROUND
1.8.0
