7 pbcoef,pwsat,pcondsat,pmpotsat,pwfc,pdg,pdzg, &
8 pdzdif,ppg,pletr,pleg,pevapcor,pf2wght, &
9 pwg,pwgi,ptg,pps,pqsat,pqsati, &
10 pdrain,phorton,pfsat,kwg_layer, &
11 kmax_layer,klayer_hort,ptopqs,pqsb,pfwtd,pwtd )
90 USE yomhook
,ONLY : lhook, dr_hook
91 USE parkind1
,ONLY : jprb
95 CHARACTER(LEN=*),
INTENT(IN) :: hrunoff
100 CHARACTER(LEN=*),
INTENT(IN) :: hhort
102 REAL,
INTENT(IN) :: ptstep
104 REAL,
DIMENSION(:),
INTENT(IN) :: pps, ppg, pletr, pleg, pevapcor, pfsat, &
119 REAL,
DIMENSION(:,:),
INTENT(IN) :: pqsat,pqsati
122 REAL,
DIMENSION(:,:),
INTENT(IN) :: ptg, pdg, pdzg, pdzdif
128 REAL,
DIMENSION(:,:),
INTENT(IN) :: pwsat, pcondsat, pf2wght, pwfc
134 REAL,
DIMENSION(:,:),
INTENT(IN) :: pbcoef,pmpotsat
138 INTEGER,
DIMENSION(:),
INTENT(IN) :: kwg_layer
141 INTEGER,
INTENT(IN) :: kmax_layer
144 INTEGER,
INTENT(IN) :: klayer_hort
146 REAL,
DIMENSION(:,:),
INTENT(INOUT) :: pwg, pwgi
150 REAL,
DIMENSION(:),
INTENT(OUT) :: pdrain, phorton
154 REAL,
DIMENSION(:,:),
INTENT(IN) :: ptopqs
157 REAL,
DIMENSION(:),
INTENT(OUT) :: pqsb
164 INTEGER :: ini, inl, idepth
166 REAL,
DIMENSION(SIZE(PDZG,1)) :: zinfiltmax, zinfiltc, zexcess, zdgn, zwgtot, zpsiwtd, zwtd
177 REAL,
DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) :: zwflux, zdfluxdt1, zdfluxdt2, zwfluxn
184 REAL,
DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) :: zpsi, zk, znu, zwsat, zhead, &
185 zvapcond, zfrz, zki, &
199 REAL,
DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) :: zamtrx, zbmtrx, zcmtrx, zfrc, zsol, &
208 REAL,
DIMENSION(SIZE(PDZG,1),SIZE(PDZG,2)) :: zinflayer
210 REAL,
PARAMETER :: zwght = 0.5
215 REAL,
PARAMETER :: zeice = 6.0
217 REAL :: zlog10, zs, zlog, zwdrain
219 REAL :: zdkdt1, zdkdt2, zdheaddt1, zdheaddt2
225 REAL(KIND=JPRB) :: zhook_handle
232 IF (lhook) CALL dr_hook(
'HYDRO_SOILDIF',0,zhook_handle)
234 ini =
SIZE(pdzg(:,:),1)
256 zcapacity(:,:) = xundef
259 zvapcond(:,:) = xundef
265 zwfluxn(:,:) = xundef
266 zdfluxdt1(:,:) = xundef
267 zdfluxdt2(:,:) = xundef
283 zwsat(jj,jl) = max(xwgmin, pwsat(jj,jl)-pwgi(jj,jl))
290 zfrz(jj,jl) = exp(zlog10*(-zeice*(pwgi(jj,jl)/(pwgi(jj,jl)+pwg(jj,jl)))))
300 IF(hrunoff==
'SGH')
THEN
301 ztopqs(:,:)=zfrz(:,:)*ptopqs(:,:)
313 IF(hrunoff/=
'WSAT'.AND.hhort/=
'SGH')
THEN
315 zinfiltmax(:) =
infmax_func(pwg,zwsat,zfrz,pcondsat,pmpotsat,pbcoef,pdzg,pdg,klayer_hort)
317 phorton(:) = (1.-pfsat(:)) * max(0.0,ppg(:)-zinfiltmax(:))
325 zinfiltc(:) = max(0.0,ppg(:)-phorton(:))*ptstep
332 zcapacity(jj,jl) = max(0.0,zwsat(jj,jl)-pwg(jj,jl))*pdzg(jj,jl)
334 zinflayer(jj,jl) = min(zinfiltc(jj),zcapacity(jj,jl))
336 pwg(jj,jl) = pwg(jj,jl)+zinflayer(jj,jl)/pdzg(jj,jl)
338 zinfiltc(jj) = zinfiltc(jj) - zinflayer(jj,jl)
352 zinfneg(jj,jl) = (min(0.0,ppg(jj))-pevapcor(jj))*pdzg(jj,jl)*pwg(jj,jl)
353 zwgtot(jj ) = zwgtot(jj)+pdzg(jj,jl)*pwg(jj,jl)
359 zinfneg(jj,jl) = zinfneg(jj,jl)/zwgtot(jj)
365 phorton(:)=(phorton(:)+zinfiltc(:)/ptstep)*xrholw
377 zs = min(1.0,pwg(jj,jl)/zwsat(jj,jl))
378 zlog = pbcoef(jj,jl)*log(zs)
379 zpsi(jj,jl) = pmpotsat(jj,jl)*exp(-zlog)
382 zlog = -zlog*(2.0+3.0/pbcoef(jj,jl))
383 zk(jj,jl) = zfrz(jj,jl)*pcondsat(jj,jl)*exp(-zlog)
394 zdgn(jj) = 0.5*(pdg(jj,idepth)+pdg(jj,idepth-1))
395 zpsiwtd(jj) = zpsi(jj,idepth)
396 IF(pwtd(jj)/=xundef)
THEN
398 zwtd(jj) = max(pdg(jj,idepth),-pwtd(jj))
400 zs = min(1.0,zwsat(jj,idepth)/pwsat(jj,idepth))
401 zlog = pbcoef(jj,idepth)*log(zs)
402 zpsiwtd(jj) = pmpotsat(jj,idepth)*exp(-zlog)
409 zvapcond(:,:) =
vapcondcf(ptg,pps,pwg,pwgi,zpsi,pwsat,pwfc,pqsat,pqsati,kwg_layer,inl)
410 zvapcond(:,:) = zfrz(:,:)*zvapcond(:,:)
423 zki(jj,jl) = sqrt(zk(jj,jl)*zk(jj,jl+1))
424 znu(jj,jl) = zki(jj,jl) + sqrt(zvapcond(jj,jl)*zvapcond(jj,jl+1))
425 zhead(jj,jl) = (zpsi(jj,jl)-zpsi(jj,jl+1))/pdzdif(jj,jl)
429 zwflux(jj,jl) = -znu(jj,jl) * zhead(jj,jl) - zki(jj,jl)
431 ELSEIF(jl==idepth)
THEN
434 zki(jj,idepth) = zk(jj,idepth)
435 znu(jj,idepth) = zk(jj,idepth) * pfwtd(jj)
436 zhead(jj,idepth) = (zpsi(jj,idepth)-zpsiwtd(jj))/(zwtd(jj)-zdgn(jj))
440 zwflux(jj,idepth) = -znu(jj,idepth) * zhead(jj,idepth) - zki(jj,idepth)
456 zdheaddt1 = -pbcoef(jj,jl )*zpsi(jj,jl )/(pwg(jj,jl )*pdzdif(jj,jl))
457 zdheaddt2 = -pbcoef(jj,jl+1)*zpsi(jj,jl+1)/(pwg(jj,jl+1)*pdzdif(jj,jl))
458 zdkdt1 = (2.*pbcoef(jj,jl )+3.)*zki(jj,jl)/(2.0*pwg(jj,jl ))
459 zdkdt2 = (2.*pbcoef(jj,jl+1)+3.)*zki(jj,jl)/(2.0*pwg(jj,jl+1))
461 zdfluxdt1(jj,jl) = -zdkdt1*zhead(jj,jl) - znu(jj,jl)*zdheaddt1 - zdkdt1
462 zdfluxdt2(jj,jl) = -zdkdt2*zhead(jj,jl) + znu(jj,jl)*zdheaddt2 - zdkdt2
463 ELSEIF(jl==idepth)
THEN
464 zdheaddt1 = -pbcoef(jj,idepth)*zpsi(jj,idepth)/(pwg(jj,idepth)*(zwtd(jj)-zdgn(jj))) &
465 +pbcoef(jj,idepth)*zpsiwtd(jj )/(zwsat(jj,idepth)*(zwtd(jj)-zdgn(jj)))
467 zdkdt1 = (2.*pbcoef(jj,idepth)+3.)*zk(jj,idepth)/pwg(jj,idepth)
470 zdfluxdt1(jj,idepth) = -zdkdt1*zhead(jj,idepth)*pfwtd(jj) - znu(jj,idepth)*zdheaddt1 - zdkdt1
471 zdfluxdt2(jj,idepth) = -zdkdt2*zhead(jj,idepth)*pfwtd(jj) + znu(jj,idepth)*zdheaddt2 - zdkdt2
480 zfrc(:,1) = zwflux(:,1) - pleg(:) - pf2wght(:,1)*pletr(:) + zinfneg(:,1) - ztopqs(:,1)
482 zbmtrx(:,1) = (pdzg(:,1)/ptstep) - zwght*zdfluxdt1(:,1)
483 zcmtrx(:,1) = -zwght*zdfluxdt2(:,1)
490 zfrc(jj,jl) = zwflux(jj,jl) - zwflux(jj,jl-1) - pf2wght(jj,jl)*pletr(jj) + zinfneg(jj,jl) - ztopqs(jj,jl)
491 zamtrx(jj,jl) = zwght*zdfluxdt1(jj,jl-1)
492 zbmtrx(jj,jl) = (pdzg(jj,jl)/ptstep) - zwght*(zdfluxdt1(jj,jl)-zdfluxdt2(jj,jl-1))
493 zcmtrx(jj,jl) = -zwght*zdfluxdt2(jj,jl)
501 CALL
tridiag_dif(zamtrx,zbmtrx,zcmtrx,zfrc,kwg_layer,inl,zsol)
514 pwg(jj,jl) = pwg(jj,jl) + zsol(jj,jl)
517 zexcess(jj) = max(0.0, pwg(jj,jl) - zwsat(jj,jl))
518 pwg(jj,jl ) = min(pwg(jj,jl),zwsat(jj,jl))
519 pwg(jj,jl+1) = pwg(jj,jl+1) + zexcess(jj)*(pdzg(jj,jl)/pdzg(jj,jl+1))
522 zwfluxn(jj,jl) = zwflux(jj,jl) + zdfluxdt1(jj,jl)*zsol(jj,jl) + zdfluxdt2(jj,jl)*zsol(jj,jl+1)
525 pqsb(jj) = pqsb(jj) + ztopqs(jj,jl)*xrholw
527 ELSEIF(jl==idepth)
THEN
530 pwg(jj,idepth) = pwg(jj,idepth) + zsol(jj,idepth)
533 zexcess(jj) = max(0.0, pwg(jj,idepth) - zwsat(jj,idepth))
534 pwg(jj,idepth) = min(pwg(jj,idepth),zwsat(jj,idepth))
535 pdrain(jj) = pdrain(jj) + zexcess(jj)*pdzg(jj,idepth)*xrholw/ptstep
538 zwfluxn(jj,idepth) = zwflux(jj,idepth) + zdfluxdt1(jj,idepth)*zsol(jj,idepth)
543 pdrain(jj) = pdrain(jj)-(zwght*zwfluxn(jj,idepth)+(1.-zwght)*zwflux(jj,idepth))*xrholw
546 pqsb(jj) = pqsb(jj) + ztopqs(jj,idepth)*xrholw
563 zexcess(jj) = max(0., xwgmin - pwg(jj,jl))
564 pwg(jj,jl) = pwg(jj,jl) + zexcess(jj)
565 pwg(jj,jl+1) = pwg(jj,jl+1) - zexcess(jj)*pdzg(jj,jl)/pdzg(jj,jl+1)
566 ELSEIF(jl==idepth.AND.pwg(jj,idepth)<xwgmin)
THEN
571 pdrain(jj) = pdrain(jj) + min(0.0,pwg(jj,idepth)-xwgmin)*pdzg(jj,idepth)*xrholw/ptstep
572 pwg(jj,idepth) = xwgmin
577 IF (lhook) CALL dr_hook(
'HYDRO_SOILDIF',1,zhook_handle)
subroutine hydro_soildif(HRUNOFF, HHORT, PTSTEP, PBCOEF, PWSAT, PCONDSAT, PMPOTSAT, PWFC, PDG, PDZG, PDZDIF, PPG, PLETR, PLEG, PEVAPCOR, PF2WGHT, PWG, PWGI, PTG, PPS, PQSAT, PQSATI, PDRAIN, PHORTON, PFSAT, KWG_LAYER, KMAX_LAYER, KLAYER_HORT, PTOPQS, PQSB, PFWTD, PWTD)