|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE HYDRO_SOILDIF(PTSTEP, & 00003 PBCOEF,PWSAT,PCONDSAT,PMPOTSAT,PWFC,PDG,PDZG,PDZDIF, & 00004 PPG,PLETR,PLEG,PEVAPCOR,PF2WGHT, & 00005 PWG,PWGI,PTG,PPS,PQSAT,PQSATI,PWDRAIN, & 00006 PDRAIN,PHORTON,HKSAT,HSOC,PWWILT,HHORT,PFSAT, & 00007 KWG_LAYER,KMAX_LAYER,KLAYER_HORT ) 00008 ! ########################################################################## 00009 ! 00010 ! 00011 !!**** *HYDRO_SOILDIF* 00012 ! 00013 !! PURPOSE 00014 !! ------- 00015 ! This subroutine solves the 1-D (z) mass conservation equation 00016 ! (mixed-form Richard's equation: tendency in terms of volumetric 00017 ! water content, gradient in terms of matric potential) 00018 ! for liquid water (using Darcy's Law for the vertical flux) 00019 ! together with the Clapp and Hornberger (1978) simplification to 00020 ! the Brooks and Corey (1966) empirical model for relating matric 00021 ! potential and hydraulic conductivity to soil water content. 00022 ! Any set of parameters can be used (eg. Clapp and Hornberger 1978; 00023 ! Cosby et al. 1984; etc.) Modifications to the equations is also 00024 ! made for the Van Genucten model/relationships. The equations 00025 ! also incorporate vapor transfer for dry soils. The soil porosity 00026 ! is modified in the presence of soil ice. Soil ice content 00027 ! is also updated herein due to changes resulting from sublimation. 00028 ! The layer averaged set of equations are time differenced 00029 ! using an implicit time scheme. 00030 ! The equations/model *assume* a heterogeneous soil texture profile, 00031 ! however, if the soil properties are homogeneous, the equations 00032 ! collapse into the standard homogeneous approach (i.e. give the 00033 ! same results as). The eqs are solved rapidly by taking advantage of the 00034 ! fact that the matrix is tridiagonal. 00035 ! Note that the exponential profile of hydraulic conductivity with soil 00036 ! depth is applied to interfacial conductivity (or interblock) 00037 ! 00038 ! 00039 !!** METHOD 00040 !! ------ 00041 ! 00042 ! Direct calculation 00043 ! 00044 !! EXTERNAL 00045 !! -------- 00046 ! 00047 ! None 00048 !! 00049 !! IMPLICIT ARGUMENTS 00050 !! ------------------ 00051 !! 00052 !! USE MODD_CST 00053 !! USE MODI_TRIDIAG_GROUND 00054 !! 00055 !! REFERENCE 00056 !! --------- 00057 !! 00058 !! Boone (2000) 00059 !! Boone et al. (2000) 00060 !! 00061 !! AUTHOR 00062 !! ------ 00063 !! A. Boone * Meteo-France * 00064 !! 00065 !! MODIFICATIONS 00066 !! ------------- 00067 !! Original 16/02/00 Boone 00068 !! Modif 04/2010 B.Decharme: geometric mean for interfacial conductivity 00069 !! Brook and Corey 00070 !! Modif 08/2011 B.Decharme: Optimization using global loops 00071 !! 10/12 B.Decharme: EVAPCOR snow correction in DIF 00072 !------------------------------------------------------------------------------- 00073 ! 00074 !* 0. DECLARATIONS 00075 ! ------------ 00076 ! 00077 USE MODD_SURF_PAR, ONLY : XUNDEF, NUNDEF 00078 USE MODD_CSTS, ONLY : XLSTT, XRHOLW, XLVTT 00079 USE MODD_ISBA_PAR, ONLY : XWGMIN 00080 ! 00081 USE MODE_HYDRO_DIF 00082 ! 00083 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00084 USE PARKIND1 ,ONLY : JPRB 00085 ! 00086 IMPLICIT NONE 00087 ! 00088 REAL, INTENT(IN) :: PTSTEP ! Model time step (s) 00089 ! 00090 REAL, DIMENSION(:), INTENT(IN) :: PPS, PPG, PLETR, PLEG, PEVAPCOR, PWDRAIN, PFSAT 00091 ! PPS = surface pressure (Pa) 00092 ! PPG = throughfall rate: 00093 ! rate at which water reaches the surface 00094 ! of the soil (from snowmelt, rain, canopy 00095 ! drip, etc...) (m/s) 00096 ! PLETR = transpiration rate (m/s) 00097 ! PLEG = bare-soil evaporation rate (m/s) 00098 ! PEVAPCOR = correction for any excess evaporation 00099 ! from snow as it completely ablates (m/s) 00100 ! PWDRAIN= minimum Wg for drainage (m3 m-3) 00101 ! PFSAT = Saturated fraction 00102 ! 00103 REAL, DIMENSION(:,:), INTENT(IN) :: PQSAT,PQSATI 00104 ! specific humidity at saturation 00105 ! 00106 REAL, DIMENSION(:,:), INTENT(IN) :: PTG, PDG, PDZG, PDZDIF 00107 ! PTG = layer-average soil temperatures (K) 00108 ! PDG = soil layer depth (m) 00109 ! PDZG= soil layer thicknesses (m) 00110 ! PDZDIF = distance between the midpoints of consecutive layers (m) 00111 ! 00112 REAL, DIMENSION(:,:), INTENT(IN) :: PWSAT, PCONDSAT, PF2WGHT, PWFC 00113 ! PWSAT = porosity profile (m3 m-3) 00114 ! PCONDSAT = hydraulic conductivity at saturation (m s-1) 00115 ! PF2WGHT = root-zone transpiration weights (-) 00116 ! PWFC = field capacity water content (m3 m-3) 00117 ! 00118 REAL, DIMENSION(:,:), INTENT(IN) :: PBCOEF,PMPOTSAT 00119 ! PMPOTSAT = matric potential at saturation (m) (BC parameters) 00120 ! PBCOEF = slope of the retention curve (-) (BC parameters) 00121 ! 00122 INTEGER, DIMENSION(:), INTENT(IN) :: KWG_LAYER 00123 ! KWG_LAYER = Number of soil moisture layers 00124 ! 00125 INTEGER, INTENT(IN) :: KMAX_LAYER 00126 ! KMAX_LAYER = Max number of soil moisture layers (DIF option) 00127 ! 00128 INTEGER, INTENT(IN) :: KLAYER_HORT! DIF optimization 00129 ! 00130 REAL, DIMENSION(:,:), INTENT(INOUT) :: PWG, PWGI 00131 ! PWG = volumetric liquid water content (m3 m-3) 00132 ! PWGI = volumetric ice content (m3 m-3) 00133 ! 00134 REAL, DIMENSION(:), INTENT(OUT) :: PDRAIN, PHORTON 00135 ! PDRAIN = drainage (flux out of model base) (kg m-2 s-1) 00136 ! PHORTON = runoff (due to saturation (lateral) (kg m-2 s-1) 00137 ! 00138 CHARACTER(LEN=*), INTENT(IN) :: HHORT ! Hortonian runoff 00139 ! 00140 CHARACTER(LEN=*), INTENT(IN) :: HKSAT ! soil hydraulic profil option 00141 ! ! 'DEF' = ISBA homogenous soil 00142 ! ! 'SGH' = ksat exponential decay 00143 ! 00144 CHARACTER(LEN=*), INTENT(IN) :: HSOC ! soil organic carbon profil option 00145 ! ! 'DEF' = ISBA homogenous soil 00146 ! ! 'SGH' = SOC profile 00147 ! 00148 REAL, DIMENSION(:,:), INTENT(IN) :: PWWILT 00149 ! PWWILT = wilting point volumetric water 00150 ! content (m3 m-3) 00151 ! 00152 ! 00153 !* 0.2 declarations of local variables 00154 ! 00155 INTEGER :: JJ, JL ! loop control 00156 ! 00157 INTEGER :: INI, INL, IDEPTH ! Number of point and grid layers 00158 ! 00159 REAL, DIMENSION(SIZE(PDZG,1)) :: ZINFILTMAX, ZINFILTC, ZEXCESS 00160 ! ZINFILTMAX = maximum allowable infiltration rate 00161 ! (from Darcy's Law) (m s-1) 00162 ! ZEXCESS = working variable: excess soil water 00163 ! which is used as a constraint 00164 ! 00165 ! 00166 REAL, DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) :: ZWFLUX, ZDFLUXDT1, ZDFLUXDT2, ZWFLUXN 00167 ! ZWFLUX = vertical soil water flux (+ up) (m s-1) 00168 ! ZDFLUXDT = total vertical flux derrivative 00169 ! ZDFLUXDT1 = vertical flux derrivative: dF_j/dw_j 00170 ! ZDFLUXDT2 = vertical flux derrivative: dF_j/dw_j+1 00171 ! ZWFLUXN = vertical soil water flux at end of time 00172 ! 00173 REAL, DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) :: ZPSI, ZK, ZNU, ZWSAT, ZHEAD, 00174 ZVAPCOND, ZFRZ, ZKI, 00175 ZCAPACITY, ZINFNEG 00176 ! ZNU = interfacial total conductivity (m s-1) 00177 ! at level z_j 00178 ! ZK = hydraulic conductivity (m s-1) 00179 ! ZHEAD = matric potential gradient (-) 00180 ! ZWSAT = ice modified soil porosity (m3 m-3) 00181 ! ZFRZ = diffusion coefficient for freezing (-) 00182 ! ZVAPCOND = vapor conductivity (m s-1) 00183 ! ZKI = interfacial hydraulic conductivity (m s-1) at level z_j 00184 ! ZCAPACITY = simple volumetric water holding capacity estimate for 00185 ! wetting front penetration (-) 00186 ! ZINFNEG = Negative infiltration (m s-1) 00187 ! 00188 REAL, DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) :: ZAMTRX, ZBMTRX, ZCMTRX, ZFRC, ZSOL, 00189 ZSGDRAIN 00190 ! ZAMTRX = leftmost diagonal element of tri-diagonal 00191 ! coefficient matrix 00192 ! ZBMTRX = center diagonal element (vector) 00193 ! ZCMTRX = rightmost diagonal element (vector) 00194 ! ZFRC = forcing function (vector) 00195 ! ZSOL = solution vector 00196 ! ZSGDRAIN = sub-grid drainge (m s-1): calibration parameter 00197 ! 00198 REAL, DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) ::ZWWILT, ZINFLAYER 00199 ! 00200 REAL, PARAMETER :: ZWGHT = 0.5 ! time scheme weight for calculating flux. 00201 ! varies from 0 (explicit time scheme) 00202 ! to 1 (backward difference implicit) 00203 ! Default is 1/2 (Crank-Nicholson) 00204 ! 00205 REAL, PARAMETER :: ZEICE = 6.0 ! Ice vertical diffusion impedence factor 00206 ! 00207 REAL :: ZLOG10, ZS, ZLOG, ZWDRAIN 00208 ! 00209 REAL :: ZWFC, ZWLIM, ZDKDT1, ZDKDT2, ZDHEADDT1, ZDHEADDT2 00210 ! ZWFC = ice modified soil field capacity (m3 m-3) 00211 ! ZDKDT1 = hydraulic conductivity derrivative w/r/t upper layer water content 00212 ! ZDKDT2 = "" lower layer water content 00213 ! ZDHEADDT1 = matric potential gradient derrivative w/r/t upper layer water content 00214 ! ZDHEADDT2 = "" lower layer water content 00215 ! 00216 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00217 ! 00218 !------------------------------------------------------------------------------- 00219 ! 0. Initialization: 00220 ! --------------- 00221 ! 00222 IF (LHOOK) CALL DR_HOOK('HYDRO_SOILDIF',0,ZHOOK_HANDLE) 00223 ! 00224 INI = SIZE(PDZG(:,:),1) 00225 INL = KMAX_LAYER 00226 ! 00227 ZLOG10 = LOG(10.0) 00228 ! 00229 PDRAIN (:) = 0.0 00230 PHORTON (:) = 0.0 00231 ZINFILTC (:) = 0.0 00232 ZEXCESS (:) = 0.0 00233 ZINFILTMAX(:) = 0.0 00234 ! 00235 ZINFNEG (:,:) = 0.0 00236 ZINFLAYER(:,:) = 0.0 00237 ZFRC (:,:) = 0.0 00238 ZFRZ (:,:) = 0.0 00239 ! 00240 ZWSAT (:,:) = XUNDEF 00241 ZCAPACITY(:,:) = XUNDEF 00242 ZPSI (:,:) = XUNDEF 00243 ZK (:,:) = XUNDEF 00244 ZVAPCOND (:,:) = XUNDEF 00245 ZKI (:,:) = XUNDEF 00246 ZSOL (:,:) = XUNDEF 00247 ZNU (:,:) = XUNDEF 00248 ZHEAD (:,:) = XUNDEF 00249 ZWFLUX (:,:) = XUNDEF 00250 ZWFLUXN (:,:) = XUNDEF 00251 ZDFLUXDT1(:,:) = XUNDEF 00252 ZDFLUXDT2(:,:) = XUNDEF 00253 ZAMTRX (:,:) = XUNDEF 00254 ZBMTRX (:,:) = XUNDEF 00255 ZCMTRX (:,:) = XUNDEF 00256 ! 00257 ZSGDRAIN (:,:) = 0.0 00258 ZWWILT (:,:) = XWGMIN 00259 IF(HKSAT=='SGH'.OR.HSOC=='SGH')THEN 00260 ZWWILT (:,:) = PWWILT(:,:) 00261 ENDIF 00262 ! 00263 DO JL=1,INL 00264 DO JJ=1,INI 00265 ! 00266 IDEPTH=KWG_LAYER(JJ) 00267 IF(JL<=IDEPTH)THEN 00268 ! 00269 ! Modification/addition of frozen soil parameters 00270 ! ----------------------------------------------- 00271 ! 00272 ! Modify soil porosity as ice assumed to become part 00273 ! of solid soil matrix (with respect to liquid flow): 00274 ZWSAT (JJ,JL) = MAX(XWGMIN, PWSAT(JJ,JL)-PWGI(JJ,JL)) 00275 ! 00276 ! Factor from (Johnsson and Lundin 1991), except here it is normalized so that it 00277 ! goes to zero in the limit as all available pore space is filled up with ice. 00278 ! For now, a simple constant is used for all soils. Further modifications 00279 ! will be made as research warrents. 00280 ! Old : 10.**(-ZEICE*(PWGI(JJ,JL)/(PWGI(JJ,JL)+PWG(JJ,JL)))) 00281 ZFRZ(JJ,JL) = EXP(ZLOG10*(-ZEICE*(PWGI(JJ,JL)/(PWGI(JJ,JL)+PWG(JJ,JL))))) 00282 ! 00283 ! Simple volumetric water holding capacity estimate for wetting front penetration 00284 ZCAPACITY(JJ,JL) = MAX(0.0,ZWSAT(JJ,JL)-PWG(JJ,JL))*PDZG(JJ,JL) 00285 00286 ! Linear (in time) sub-grid drainage term (input in m s-1 herein) 00287 ! --------------------------------------------------------------- 00288 ! 00289 ! (very simple way to take into account sub-grid water table effects, pores, etc...): 00290 ! Can be calibrated. Also, make this effect vanish/cease for 00291 ! extremely dry soil layers. This term/option is OFF if WDRAIN = 0. 00292 ! 00293 IF(PWDRAIN(JJ) > 0.0)THEN 00294 ZWFC = PWFC (JJ,JL)*ZWSAT(JJ,JL)/PWSAT(JJ,JL) 00295 ZWLIM = MAX(XWGMIN,ZWWILT(JJ,JL)*ZWSAT(JJ,JL)/PWSAT(JJ,JL)) 00296 ZS = MIN(1.0,(ZWFC+PWDRAIN(JJ))/ZWSAT(JJ,JL)) 00297 ZLOG = (2.*PBCOEF(JJ,JL)+3.0)*LOG(ZS) 00298 ZWDRAIN = ZFRZ(JJ,JL)*PCONDSAT(JJ,JL) * EXP(ZLOG) 00299 ZSGDRAIN(JJ,JL) = ZWDRAIN * MAX(0.0, MIN(ZWFC,PWG(JJ,JL))-ZWLIM)/(ZWFC-ZWLIM) 00300 ENDIF 00301 ! 00302 ENDIF 00303 ! 00304 ENDDO 00305 ENDDO 00306 ! 00307 ! 1. Infiltration at "t" 00308 ! ------------------- 00309 ! 00310 ! Surface flux term (m s-1): Infiltration (and surface runoff) 00311 ! Surface fluxes are limited a Green-Ampt approximation from Abramopoulos et al 00312 ! (1988) and Entekhabi and Eagleson (1989). 00313 ! Note : when Horton option is used, infiltration already calculated in hydro_sgh 00314 IF(HHORT/='SGH')THEN 00315 ! Green-Ampt approximation for maximum infiltration (derived form) 00316 ZINFILTMAX(:) = INFMAX_FUNC(PWG,ZWSAT,ZFRZ,PCONDSAT,PMPOTSAT,PBCOEF,PDZG,PDG,KLAYER_HORT) 00317 ! Fast(temporal)-response runoff (surface excess) (m s-1): 00318 PHORTON (:) = (1.-PFSAT(:)) * MAX(0.0,PPG(:)-ZINFILTMAX(:)) 00319 ENDIF 00320 ! 00321 ! 00322 ! 2. Initialise soil moisture profile according to the sink linear terms at "t" 00323 ! ---------------------------------------------------------------------------- 00324 ! 00325 !Surface cumulative infiltration (m) 00326 ZINFILTC(:) = MAX(0.0,PPG(:)-PHORTON(:))*PTSTEP 00327 ! 00328 DO JL=1,INL 00329 DO JJ=1,INI 00330 IDEPTH=KWG_LAYER(JJ) 00331 IF(JL<=IDEPTH)THEN 00332 ! Infiltration terms (m) : 00333 ZINFLAYER(JJ,JL) = MIN(ZINFILTC(JJ),ZCAPACITY(JJ,JL)) 00334 ! Soil moisture (m3/m3) : 00335 PWG (JJ,JL) = PWG(JJ,JL)+ZINFLAYER(JJ,JL)/PDZG(JJ,JL) 00336 ! Put remainding infiltration into the next layer (m) 00337 ZINFILTC(JJ) = ZINFILTC(JJ) - ZINFLAYER(JJ,JL) 00338 ! Possible negative infiltration (m s-1) 00339 ZINFNEG(JJ,JL) = (MIN(0.0,PPG(JJ))+PEVAPCOR(JJ))*PDZG(JJ,JL)/PDG(JJ,IDEPTH) 00340 ENDIF 00341 ENDDO 00342 ENDDO 00343 ! 00344 !Fast(temporal)-response runoff (surface excess) (kg m2 s-1): 00345 !special case : if infiltration > total soil capacity 00346 PHORTON(:)=(PHORTON(:)+ZINFILTC(:)/PTSTEP)*XRHOLW 00347 ! 00348 ! 00349 ! 3. Initialise matric potential and hydraulic conductivity at "t" 00350 ! ------------------------------------------------------------- 00351 ! 00352 DO JL=1,INL 00353 DO JJ=1,INI 00354 IDEPTH=KWG_LAYER(JJ) 00355 IF(JL<=IDEPTH)THEN 00356 ! Matric potential (m) : 00357 ! psi=mpotsat*(w/wsat)**(-bcoef) 00358 ZS = MIN(1.0,PWG(JJ,JL)/ZWSAT(JJ,JL)) 00359 ZLOG = PBCOEF(JJ,JL)*LOG(ZS) 00360 ZPSI(JJ,JL) = PMPOTSAT(JJ,JL)*EXP(-ZLOG) 00361 ! Hydraulic conductivity from matric potential (m s-1): 00362 ! k=frz*condsat*(psi/mpotsat)**(-2-3/bcoef) 00363 ZLOG = -ZLOG*(2.0+3.0/PBCOEF(JJ,JL)) 00364 ZK(JJ,JL) = ZFRZ(JJ,JL)*PCONDSAT(JJ,JL)*EXP(-ZLOG) 00365 ENDIF 00366 ENDDO 00367 ENDDO 00368 ! 00369 ! 00370 ! 4. Vapor diffusion conductivity (m s-1) 00371 ! ------------------------------------ 00372 ! 00373 ZVAPCOND(:,:) = VAPCONDCF(PTG,PPS,PWG,PWGI,ZPSI,PWSAT,PWFC,PQSAT,PQSATI,KWG_LAYER,INL) 00374 ZVAPCOND(:,:) = ZFRZ(:,:)*ZVAPCOND(:,:) 00375 ! 00376 ! 5. Linearized water flux: values at "t" 00377 ! ------------------------------------ 00378 ! calculate flux at the beginning of the time step: 00379 ! 00380 DO JL=1,INL 00381 DO JJ=1,INI 00382 ! 00383 IDEPTH=KWG_LAYER(JJ) 00384 IF(JL<IDEPTH)THEN 00385 ! 00386 ! Total interfacial conductivity (m s-1) And Potential gradient (dimensionless): 00387 ZKI (JJ,JL) = SQRT(ZK(JJ,JL)*ZK(JJ,JL+1)) 00388 ZNU (JJ,JL) = ZKI(JJ,JL) + SQRT(ZVAPCOND(JJ,JL)*ZVAPCOND(JJ,JL+1)) 00389 ZHEAD(JJ,JL) = (ZPSI(JJ,JL)-ZPSI(JJ,JL+1))/PDZDIF(JJ,JL) 00390 ! 00391 ! Total Sub-surface soil water fluxes (m s-1): (+ up, - down) using Darcy's 00392 ! Law with an added linear drainage term: 00393 ZWFLUX(JJ,JL) = -ZNU(JJ,JL) *ZHEAD(JJ,JL) - ZKI(JJ,JL) - ZSGDRAIN(JJ,JL) 00394 ! 00395 ELSEIF(JL==IDEPTH)THEN !Last layers 00396 ! 00397 ! Total interfacial conductivity (m s-1) And Potential gradient (dimensionless): 00398 ZKI (JJ,IDEPTH) = ZK(JJ,IDEPTH) 00399 ZNU (JJ,IDEPTH) = 0.0 00400 ZHEAD(JJ,IDEPTH) = 0.0 00401 ! 00402 ! Total Sub-surface soil water fluxes (m s-1): (+ up, - down) using Darcy's 00403 ! Law with an added linear drainage term: 00404 ZWFLUX(JJ,IDEPTH) = -ZNU(JJ,IDEPTH) *ZHEAD(JJ,IDEPTH) - ZKI(JJ,IDEPTH) - ZSGDRAIN(JJ,IDEPTH) 00405 ! 00406 ENDIF 00407 ! 00408 ENDDO 00409 ENDDO 00410 ! 00411 ! 00412 ! 7. Linearized water flux: values at "t+dt" 00413 ! --------------------------------------- 00414 ! Flux Derrivative terms, see A. Boone thesis (Annexe E). 00415 ! 00416 DO JL=1,INL 00417 DO JJ=1,INI 00418 IDEPTH=KWG_LAYER(JJ) 00419 IF(JL<IDEPTH)THEN 00420 ZDHEADDT1 = -PBCOEF(JJ,JL )*ZPSI(JJ,JL )/(PWG(JJ,JL )*PDZDIF(JJ,JL)) 00421 ZDHEADDT2 = -PBCOEF(JJ,JL+1)*ZPSI(JJ,JL+1)/(PWG(JJ,JL+1)*PDZDIF(JJ,JL)) 00422 ZDKDT1 = (2.*PBCOEF(JJ,JL )+3.)*ZKI(JJ,JL)/(2.0*PWG(JJ,JL )) 00423 ZDKDT2 = (2.*PBCOEF(JJ,JL+1)+3.)*ZKI(JJ,JL)/(2.0*PWG(JJ,JL+1)) 00424 ! Total Flux derrivative terms: 00425 ZDFLUXDT1(JJ,JL) = -ZDKDT1*ZHEAD(JJ,JL) - ZNU(JJ,JL)*ZDHEADDT1 - ZDKDT1 00426 ZDFLUXDT2(JJ,JL) = -ZDKDT2*ZHEAD(JJ,JL) + ZNU(JJ,JL)*ZDHEADDT2 - ZDKDT2 00427 ELSEIF(JL==IDEPTH)THEN !Last layers 00428 ZDHEADDT1 = 0.0 00429 ZDHEADDT2 = 0.0 00430 ZDKDT1 = (2.*PBCOEF(JJ,IDEPTH)+3.)*ZK(JJ,IDEPTH)/PWG(JJ,IDEPTH) 00431 ZDKDT2 = 0.0 00432 ! Total Flux derrivative terms: 00433 ZDFLUXDT1(JJ,IDEPTH) = -ZDKDT1*ZHEAD(JJ,IDEPTH) - ZNU(JJ,IDEPTH)*ZDHEADDT1 - ZDKDT1 00434 ZDFLUXDT2(JJ,IDEPTH) = -ZDKDT2*ZHEAD(JJ,IDEPTH) + ZNU(JJ,IDEPTH)*ZDHEADDT2 - ZDKDT2 00435 ENDIF 00436 ENDDO 00437 ENDDO 00438 ! 00439 ! 8. Jacobian Matrix coefficients and Forcing function 00440 ! ------------------------------------------------- 00441 ! 00442 !surface layer: 00443 ZFRC (:,1) = ZWFLUX(:,1) - PLEG(:) - PF2WGHT(:,1)*PLETR(:) + ZINFNEG(:,1) 00444 ZAMTRX(:,1) = 0.0 00445 ZBMTRX(:,1) = (PDZG(:,1)/PTSTEP) - ZWGHT*ZDFLUXDT1(:,1) 00446 ZCMTRX(:,1) = -ZWGHT*ZDFLUXDT2(:,1) 00447 ! 00448 !Other sub-surface layers: 00449 DO JL=2,INL 00450 DO JJ=1,INI 00451 IDEPTH=KWG_LAYER(JJ) 00452 IF(JL<=IDEPTH)THEN 00453 ZFRC (JJ,JL) = ZWFLUX (JJ,JL) - ZWFLUX(JJ,JL-1) - PF2WGHT(JJ,JL)*PLETR(JJ) + ZINFNEG(JJ,JL) 00454 ZAMTRX(JJ,JL) = ZWGHT*ZDFLUXDT1(JJ,JL-1) 00455 ZBMTRX(JJ,JL) = (PDZG(JJ,JL)/PTSTEP) - ZWGHT*(ZDFLUXDT1(JJ,JL)-ZDFLUXDT2(JJ,JL-1)) 00456 ZCMTRX(JJ,JL) = -ZWGHT*ZDFLUXDT2(JJ,JL) 00457 ENDIF 00458 ENDDO 00459 ENDDO 00460 ! 00461 ! Solve Matrix Equation: tridiagonal system: solve for soil 00462 ! water (volumetric water content) tendencies: 00463 ! 00464 CALL TRIDIAG_DIF(ZAMTRX,ZBMTRX,ZCMTRX,ZFRC,KWG_LAYER,INL,ZSOL) 00465 ! 00466 ! 9. Final calculations and diagnostics: 00467 ! ----------------------------------- 00468 ! 00469 ! 00470 DO JL=1,INL 00471 DO JJ=1,INI 00472 ! 00473 IDEPTH=KWG_LAYER(JJ) 00474 IF(JL<IDEPTH)THEN 00475 ! 00476 ! Update liquid water content (m3 m-3): 00477 PWG(JJ,JL) = PWG(JJ,JL) + ZSOL(JJ,JL) 00478 ! 00479 ! Supersaturated drainage (kg m-2 s-1): 00480 ZEXCESS(JJ) = MAX(0.0, PWG(JJ,JL) - ZWSAT(JJ,JL)) 00481 PWG(JJ,JL ) = MIN(PWG(JJ,JL),ZWSAT(JJ,JL)) 00482 PWG(JJ,JL+1) = PWG(JJ,JL+1) + ZEXCESS(JJ)*(PDZG(JJ,JL)/PDZG(JJ,JL+1)) 00483 ! 00484 ! final fluxes (at end of time step) (m s-1): 00485 ZWFLUXN(JJ,JL) = ZWFLUX(JJ,JL) + ZDFLUXDT1(JJ,JL)*ZSOL(JJ,JL) + ZDFLUXDT2(JJ,JL)*ZSOL(JJ,JL+1) 00486 ! 00487 ELSEIF(JL==IDEPTH)THEN 00488 ! 00489 ! Update liquid water content (m3 m-3): 00490 PWG(JJ,IDEPTH) = PWG(JJ,IDEPTH) + ZSOL(JJ,IDEPTH) 00491 ! 00492 ! Supersaturated drainage (kg m-2 s-1): 00493 ZEXCESS(JJ) = MAX(0.0, PWG(JJ,IDEPTH) - ZWSAT(JJ,IDEPTH)) 00494 PWG(JJ,IDEPTH) = MIN(PWG(JJ,IDEPTH),ZWSAT(JJ,IDEPTH)) 00495 PDRAIN (JJ) = PDRAIN(JJ) + ZEXCESS(JJ)*PDZG(JJ,IDEPTH)*XRHOLW/PTSTEP 00496 ! 00497 ! final fluxes (at end of time step) (m s-1): 00498 ZWFLUXN(JJ,IDEPTH) = ZWFLUX(JJ,IDEPTH) + ZDFLUXDT1(JJ,IDEPTH)*ZSOL(JJ,IDEPTH) 00499 ! 00500 ! Drainage or baseflow out of bottom of model (slow time response) (kg m-2 s-1): 00501 ! Final fluxes (if needed) over the time step (kg m-2 s-1) 00502 ! would be calculated as (for water budget checks) as F = [ wgt*F(t+dt) + (1.-wgt)*F(t) ]*XRHOLW 00503 PDRAIN (JJ) = PDRAIN(JJ)-(ZWGHT*ZWFLUXN(JJ,IDEPTH)+(1.-ZWGHT)*ZWFLUX(JJ,IDEPTH))*XRHOLW 00504 ! 00505 ENDIF 00506 ! 00507 ENDDO 00508 ENDDO 00509 ! 00510 ! Possible correction/Constraint application: 00511 ! 00512 DO JL=1,INL 00513 DO JJ=1,INI 00514 IDEPTH=KWG_LAYER(JJ) 00515 IF(JL<IDEPTH)THEN 00516 ! if the soil water happens to fall below the minimum, then 00517 ! extract needed water from the layer below: this should 00518 ! generally be non-existant: but added to ensure conservation 00519 ! even for the most extreme events. 00520 ZEXCESS(JJ) = MAX(0., XWGMIN - PWG(JJ,JL)) 00521 PWG(JJ,JL) = PWG(JJ,JL) + ZEXCESS(JJ) 00522 PWG(JJ,JL+1) = PWG(JJ,JL+1) - ZEXCESS(JJ)*PDZG(JJ,JL)/PDZG(JJ,JL+1) 00523 ELSEIF(JL==IDEPTH.AND.PWG(JJ,IDEPTH)<XWGMIN)THEN 00524 ! NOTE, negative moisture can arise for *completely* dry/dessicated soils 00525 ! owing to the above check because vertical fluxes 00526 ! can be *very* small but nonzero. Here correct owing to 00527 ! small numerical drainage. 00528 PDRAIN(JJ) = PDRAIN(JJ) + MIN(0.0,PWG(JJ,IDEPTH)-XWGMIN)*PDZG(JJ,IDEPTH)*XRHOLW/PTSTEP 00529 PWG(JJ,IDEPTH) = XWGMIN 00530 ENDIF 00531 ENDDO 00532 ENDDO 00533 ! 00534 IF (LHOOK) CALL DR_HOOK('HYDRO_SOILDIF',1,ZHOOK_HANDLE) 00535 ! 00536 END SUBROUTINE HYDRO_SOILDIF 00537 00538 00539 00540 00541 00542
1.8.0
