|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE TRIDIAG_SURF(PVARM,PF,PDFDDTDZ,PEXT,PDEXTDV,PTSTEP, & 00003 PDZZ,PDZM,PVARP,OIMPL,PALFA,PBETA ) 00004 ! ################################################# 00005 ! 00006 ! 00007 !!**** *TRIDIAG_SURF* - routine to solve a time implicit scheme 00008 !! 00009 !! 00010 !! PURPOSE 00011 !! ------- 00012 ! The purpose of this routine is to give a field PVARP at t+1, by 00013 ! solving an implicit TRIDIAGonal system obtained by the 00014 ! discretization of the vertical turbulent diffusion. It should be noted 00015 ! that the degree of implicitness can be varied (ZIMPL parameter) and that 00016 ! the function of F(dT/dz) must have been linearized. 00017 ! PVARP is localized at a mass point. 00018 ! 00019 !!** METHOD 00020 !! ------ 00021 !! 00022 !! [T(+) - T(-)]/Dt = -d{ F + dF/d(dT/dz) * impl * [dT/dz(+) - dT/dz(-)] }/dz 00023 !! + Ext + dExt/dT * impl [ T(+) - T(-)] 00024 !! 00025 !! It is discretized as follows: 00026 !! 00027 !! PDZM(k)*PVARP(k)/PTSTEP 00028 !! = 00029 !! PDZM(k)*PVARM(k)/PTSTEP 00030 !! - (PDZM(k+1)+PDZM(k) )/2. * PF(k+1)/PDZZ(k+1) 00031 !! + (PDZM(k) +PDZM(k-1))/2. * PF(k) /PDZZ(k) 00032 !! + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2 00033 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARP(k+1)/PDZZ(k+1)**2 00034 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2 00035 !! + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARP(k) /PDZZ(k+1)**2 00036 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2 00037 !! + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARP(k) /PDZZ(k)**2 00038 !! + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2 00039 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARP(k-1)/PDZZ(k)**2 00040 !! + PDZM(k) * PEXT(k) 00041 !! + PDZM(k) * ZIMPL* PDEXTDV(k) * PVARP(k) 00042 !! - PDZM(k) * ZIMPL* PDEXTDV(k) * PVARM(k) 00043 !! 00044 !! 00045 !! 00046 !! The system to solve is: 00047 !! 00048 !! A*PVARP(k-1) + B*PVARP(k) + C*PVARP(k+1) = Y(k) 00049 !! 00050 !! 00051 !! The RHS of the linear system in PVARP writes: 00052 !! 00053 !! y(k) = PDZM(k)*PVARM(k)/PTSTEP 00054 !! - (PDZM(k+1)+PDZM(k) )/2. * PF(k+1)/PDZZ(k+1) 00055 !! + (PDZM(k) +PDZM(k-1))/2. * PF(k) /PDZZ(k) 00056 !! + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2 00057 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2 00058 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2 00059 !! + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2 00060 !! + PDZM(k) * PEXT(k) 00061 !! - PDZM(k) * PDEXTDV(k) * PVARM(k) 00062 !! 00063 !! 00064 !! Then, the classical TRIDIAGonal algorithm is used to invert the 00065 !! implicit operator. Its matrix is given by: 00066 !! 00067 !! ( b(ikb) c(ikb) 0 0 0 0 0 0 ) 00068 !! ( 0 a(ikb+1) b(ikb+1) c(ikb+1) 0 ... 0 0 0 ) 00069 !! ( 0 0 a(ikb+2) b(ikb+2) c(ikb+2). 0 0 0 ) 00070 !! ....................................................................... 00071 !! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 ) 00072 !! ....................................................................... 00073 !! ( 0 0 0 0 0 ...a(ike-1) b(ike-1) c(ike-1)) 00074 !! ( 0 0 0 0 0 ... 0 a(ike) b(ike) ) 00075 !! 00076 !! ikb and ike represent the first and the last inner mass levels of the 00077 !! model. The coefficients are: 00078 !! 00079 !! a(k) = + (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) /PDZZ(k)**2 00080 !! b(k) = PDZM(k) / PTSTEP 00081 !! - (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2 00082 !! - (PDZM(k) +PDZM(k-1))/2. * ZIMPL* PDFDDTDZ(k) /PDZZ(k)**2 00083 !! - PDZM(k) * PDEXTDV(k) 00084 !! c(k) = + (PDZM(k+1)+PDZM(k) )/2. * ZIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2 00085 !! 00086 !! for all k /= ikb or ike 00087 !! 00088 !! 00089 !! b(ikb) = PDZM(ikb) / PTSTEP 00090 !! -(PDZM(ikb+1)+PDZM(ikb))/2.*ZIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2 00091 !! c(ikb) = +(PDZM(ikb+1)+PDZM(ikb))/2.*ZIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2 00092 !! 00093 !! b(ike) = PDZM(ike) / PTSTEP 00094 !! -(PDZM(ike)+PDZM(ike-1))/2.*ZIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2 00095 !! a(ike) = +(PDZM(ike)+PDZM(ike-1))/2.*ZIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2 00096 !! 00097 !! 00098 !! EXTERNAL 00099 !! -------- 00100 !! 00101 !! NONE 00102 !! 00103 !! IMPLICIT ARGUMENTS 00104 !! ------------------ 00105 !! 00106 !! REFERENCE 00107 !! --------- 00108 !! Press et al: Numerical recipes (1986) Cambridge Univ. Press 00109 !! 00110 !! AUTHOR 00111 !! ------ 00112 !! V. Masson * Meteo-France * 00113 !! 00114 !! MODIFICATIONS 00115 !! ------------- 00116 !! Original 04/2003 (from tridiag.f90) 00117 !! 09/2007 implicitation with surface flux at lowest level 00118 !! --------------------------------------------------------------------- 00119 ! 00120 !* 0. DECLARATIONS 00121 ! 00122 ! 00123 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00124 USE PARKIND1 ,ONLY : JPRB 00125 ! 00126 IMPLICIT NONE 00127 ! 00128 ! 00129 !* 0.1 declarations of arguments 00130 ! 00131 REAL, DIMENSION(:,:), INTENT(IN) :: PVARM ! variable at t-1 at mass point 00132 REAL, DIMENSION(:,:), INTENT(IN) :: PF ! flux in dV/dt=-dF/dz at flux point 00133 REAL, DIMENSION(:,:), INTENT(IN) :: PDFDDTDZ! dF/d(dV/dz) at flux point 00134 REAL, DIMENSION(:,:), INTENT(IN) :: PEXT ! External forces in dT/dt=Ext 00135 ! ! (Except surf. flux) at mass point 00136 REAL, DIMENSION(:,:), INTENT(IN) :: PDEXTDV ! dExt/dV 00137 ! ! at mass point 00138 REAL, INTENT(IN) :: PTSTEP ! time step 00139 REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ ! Dz at flux point 00140 REAL, DIMENSION(:,:), INTENT(IN) :: PDZM ! Dz at mass point 00141 ! 00142 REAL, DIMENSION(:,:), INTENT(OUT):: PVARP ! variable at t+1 at mass point 00143 LOGICAL, OPTIONAL, INTENT(IN) :: OIMPL ! true if implicit coupling 00144 ! ! information between lvl 1 and 00145 ! ! surface flux scheme must be 00146 ! ! returned. 00147 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PALFA ! V+(1) = alfa F(1) + beta 00148 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PBETA ! V+(1) = alfa F(1) + beta 00149 ! 00150 ! 00151 !* 0.2 declarations of local variables 00152 ! 00153 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZVARP 00154 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZDZ_DFDDTDZ_O_DZ2 00155 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZA, ZB, ZC 00156 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2)) :: ZY ,ZGAM 00157 ! RHS of the equation, 3D work array 00158 REAL, DIMENSION(SIZE(PVARM,1)) :: ZBET 00159 ! 2D work array 00160 INTEGER :: JK ! loop counter 00161 INTEGER :: IK ! vertical limits 00162 ! 00163 LOGICAL :: GIMPL 00164 REAL :: ZIMPL ! implicitation coefficient 00165 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00166 ! ! for SBL scheme solving 00167 ! --------------------------------------------------------------------------- 00168 ! 00169 !* 1. Preliminaries 00170 ! ------------- 00171 ! 00172 IF (LHOOK) CALL DR_HOOK('TRIDIAG_SURF',0,ZHOOK_HANDLE) 00173 ZIMPL = 1 00174 ! 00175 IK = SIZE(PF,2) 00176 ! 00177 ZDZ_DFDDTDZ_O_DZ2 = PDFDDTDZ/PDZZ 00178 ! 00179 ZA=0. 00180 ZB=0. 00181 ZC=0. 00182 ZY=0. 00183 ! 00184 IF (PRESENT(OIMPL)) THEN 00185 GIMPL = OIMPL 00186 ELSE 00187 GIMPL = .FALSE. 00188 END IF 00189 ! 00190 !* 2. COMPUTE THE RIGHT HAND SIDE 00191 ! --------------------------- 00192 ! 00193 ZY(:,1) = PDZM(:,1)*PVARM(:,1)/PTSTEP & 00194 - PF(:,2) & 00195 + PF(:,1) & 00196 - ZDZ_DFDDTDZ_O_DZ2(:,1) * ZIMPL * PVARM(:,1)& 00197 + ZDZ_DFDDTDZ_O_DZ2(:,2) * ZIMPL * PVARM(:,2)& 00198 - ZDZ_DFDDTDZ_O_DZ2(:,2) * ZIMPL * PVARM(:,1)& 00199 + PDZM(:,1)*PEXT(:,1) & 00200 - PDZM(:,1)*PDEXTDV(:,1) * ZIMPL * PVARM(:,1) 00201 ! 00202 DO JK=2,IK-1 00203 ZY(:,JK) = PDZM(:,JK)*PVARM(:,JK)/PTSTEP & 00204 - PF(:,JK+1) & 00205 + PF(:,JK ) & 00206 + ZDZ_DFDDTDZ_O_DZ2(:,JK+1) * ZIMPL * PVARM(:,JK+1) & 00207 - ZDZ_DFDDTDZ_O_DZ2(:,JK+1) * ZIMPL * PVARM(:,JK ) & 00208 - ZDZ_DFDDTDZ_O_DZ2(:,JK ) * ZIMPL * PVARM(:,JK ) & 00209 + ZDZ_DFDDTDZ_O_DZ2(:,JK ) * ZIMPL * PVARM(:,JK-1) & 00210 + PDZM(:,JK)*PEXT(:,JK) & 00211 - PDZM(:,JK)*PDEXTDV(:,JK) * ZIMPL * PVARM(:,JK) 00212 END DO 00213 ! 00214 !* upper level is fixed (atmospheric forcing) 00215 ! turbulent flux divergence is supposed equal to ZERO (inertial sublayer). 00216 ! other terms are kept 00217 ! 00218 ZY(:,IK) = PDZM(:,IK)*PVARM(:,IK)/PTSTEP & 00219 + PDZM(:,IK)*PEXT(:,IK) & 00220 - PDZM(:,IK)*PDEXTDV(:,IK) * ZIMPL * PVARM(:,IK) 00221 ! 00222 ! 00223 !* 3. INVERSION OF THE TRIDIAGONAL SYSTEM 00224 ! ----------------------------------- 00225 ! 00226 ! 00227 !* 3.1 arrays A, B, C 00228 ! -------------- 00229 ! 00230 ZB(:,1) = PDZM(:,1)/PTSTEP & 00231 - ZDZ_DFDDTDZ_O_DZ2(:,1) * ZIMPL & 00232 - ZDZ_DFDDTDZ_O_DZ2(:,2) * ZIMPL & 00233 - PDZM(:,1)*PDEXTDV(:,1) 00234 ZC(:,1) = ZDZ_DFDDTDZ_O_DZ2(:,2) * ZIMPL 00235 00236 DO JK=2,IK-1 00237 ZA(:,JK) = ZDZ_DFDDTDZ_O_DZ2(:,JK ) * ZIMPL 00238 ZB(:,JK) = PDZM(:,JK)/PTSTEP & 00239 - ZDZ_DFDDTDZ_O_DZ2(:,JK+1) * ZIMPL & 00240 - ZDZ_DFDDTDZ_O_DZ2(:,JK ) * ZIMPL & 00241 - PDZM(:,JK)*PDEXTDV(:,JK) 00242 ZC(:,JK) = ZDZ_DFDDTDZ_O_DZ2(:,JK+1) * ZIMPL 00243 END DO 00244 00245 !* upper level is fixed (atmospheric forcing) 00246 ! turbulent flux divergence is supposed equal to ZERO (inertial sublayer). 00247 ! other terms are kept 00248 ZA(:,IK) = 0. 00249 ZB(:,IK) = PDZM(:,IK)/PTSTEP & 00250 - PDZM(:,IK)*PDEXTDV(:,IK) 00251 ! 00252 !* 3.2 going down 00253 ! ---------- 00254 ! 00255 ZBET(:) = ZB(:,IK) ! bet = b(ik) 00256 ZVARP(:,IK) = ZY(:,IK) / ZBET(:) 00257 00258 ! 00259 DO JK = IK-1,2,-1 00260 ZGAM(:,JK) = ZA(:,JK+1) / ZBET(:) 00261 ! gam(k) = c(k-1) / bet 00262 ZBET(:) = ZB(:,JK) - ZC(:,JK) * ZGAM(:,JK) 00263 ! bet = b(k) - a(k)* gam(k) 00264 ZVARP(:,JK)= ( ZY(:,JK) - ZC(:,JK) * ZVARP(:,JK+1) ) / ZBET(:) 00265 ! res(k) = (y(k) -a(k)*res(k-1))/ bet 00266 END DO 00267 ! special treatment for the lowest level 00268 ZGAM(:,1) = ZA(:,2) / ZBET(:) 00269 ! gam(k) = c(k-1) / bet 00270 ZBET(:) = ZB(:,1) - ZC(:,1) * ZGAM(:,1) 00271 ! bet = b(k) - a(k)* gam(k) 00272 ZVARP(:,1)= ( ZY(:,1) - ZC(:,1) * ZVARP(:,2) ) / ZBET(:) 00273 ! res(k) = (y(k) -a(k)*res(k-1))/ bet 00274 ! 00275 ! 00276 !* 3.3 going up 00277 ! -------- 00278 ! 00279 DO JK = 2,IK 00280 ZVARP(:,JK) = ZVARP(:,JK) - ZGAM(:,JK-1) * ZVARP(:,JK-1) 00281 END DO 00282 00283 ! 00284 ! 00285 !* 3.4 updates prognostic variable 00286 ! --------------------------- 00287 ! 00288 IF (.NOT. GIMPL) THEN 00289 PVARP=ZVARP 00290 ELSE 00291 PALFA=1./ZBET(:) 00292 PBETA=ZVARP(:,1) 00293 PVARP=PVARM 00294 END IF 00295 IF (LHOOK) CALL DR_HOOK('TRIDIAG_SURF',1,ZHOOK_HANDLE) 00296 ! 00297 !------------------------------------------------------------------------------- 00298 ! 00299 END SUBROUTINE TRIDIAG_SURF
1.8.0
