6 SUBROUTINE soil( HC1DRY, HSCOND, HSNOW_ISBA, OGLACIER, &
9 pc1sat, pc2ref, pacoef, ppcoef, pcv, &
10 ppsn, ppsng, ppsnv, pffg, pffv, pff, &
11 pcg, pc1, pc2, pwgeq, pct, pcs, pfrozen1, &
13 phcapsoilz, pconddryz, pcondsldz, &
14 pbcoef, pwsat, pwwilt, &
15 hksat, pcondsat, pffg_nosnow, pffv_nosnow )
77 USE modd_csts, ONLY : xpi, xci, xrholi, xday, xcl, xrholw, xcondi
83 USE yomhook
,ONLY : lhook, dr_hook
84 USE parkind1
,ONLY : jprb
91 CHARACTER(LEN=*),
INTENT(IN) :: hc1dry
96 CHARACTER(LEN=*),
INTENT(IN) :: hscond
103 CHARACTER(LEN=*),
INTENT(IN) :: hsnow_isba
108 LOGICAL,
INTENT(IN) :: oglacier
113 REAL,
DIMENSION(:),
INTENT(IN) :: psnowrhom
117 REAL,
INTENT(IN) :: pcgmax
120 REAL,
DIMENSION(:),
INTENT(IN) :: pveg, pcgsat, pc1sat
121 REAL,
DIMENSION(:),
INTENT(IN) :: pc2ref, pacoef, ppcoef, pcv
131 REAL,
DIMENSION(:),
INTENT(IN) :: ppsn, ppsng, ppsnv
138 REAL,
DIMENSION(:),
INTENT(IN) :: phcapsoilz, pconddryz, pcondsldz
145 REAL,
DIMENSION(:),
INTENT(IN) :: pbcoef, pwsat, pwwilt
150 REAL,
DIMENSION(:),
INTENT(IN) :: ptg
153 REAL,
DIMENSION(:,:),
INTENT(IN) :: pwg, pwgi
157 REAL,
DIMENSION(:),
INTENT(OUT) :: pcg, pc1, pc2, pwgeq, pct
158 REAL,
DIMENSION(:),
INTENT(OUT) :: pcs, pfrozen1
168 CHARACTER(LEN=*),
INTENT(IN) :: hksat
172 REAL,
DIMENSION(:,:),
INTENT(IN) :: pcondsat
176 REAL,
DIMENSION(:),
INTENT(IN) :: pffv, pffg, pff, pffg_nosnow, pffv_nosnow
184 REAL,
DIMENSION(SIZE(PVEG)) :: zlams, &
187 zcw1max, zx2, zy1, zy2, &
188 zlymy1, zza, zzb, zdelta, &
212 REAL,
DIMENSION(SIZE(PVEG)) :: zfrozen2, zunfrozen2, zcondsat, zsatdeg, zkersten, &
223 REAL,
DIMENSION(SIZE(PVEG)) :: zwg2
226 REAL,
DIMENSION(SIZE(PVEG)) :: zcf
227 REAL,
DIMENSION(SIZE(PVEG)) :: zff
230 REAL(KIND=JPRB) :: zhook_handle
233 IF (lhook) CALL dr_hook(
'SOIL',0,zhook_handle)
267 WHERE (pwgi(:,1) + pwg(:,1) .NE. 0.)
268 pfrozen1(:) = pwgi(:,1) / (pwgi(:,1) + pwg(:,1))
273 zwsat(jj) = max(pwsat(jj) - pwgi(jj,2),xwgmin)
275 zwsat1(jj) = max(pwsat(jj) - pwgi(jj,1),xwgmin)
277 zwwilt(jj) = pwwilt(jj) * (zwsat1(jj) / pwsat(jj))
285 IF(hscond ==
'NP89')
THEN
301 pcg(:) = (1.-pwgi(:,2)) * pcgsat(:) * ( zwsat(:)/pwg(:,2) ) &
302 **( 0.5*pbcoef(:)/log(10.) ) &
303 + pwgi(:,2) * 2. * sqrt(xpi/(xcondi*xci*xrholi*xday))
320 zfrozen2(jj) = pwgi(jj,2)/(pwgi(jj,2) + pwg(jj,2))
324 zunfrozen2(jj) = (1.0-zfrozen2(jj))*pwsat(jj)
328 zcondsat(jj) = (pcondsldz(jj)**(1.0-pwsat(jj)))* &
329 (xcondi**(pwsat(jj)-zunfrozen2(jj)))* &
330 (xcondwtr**zunfrozen2(jj))
334 zsatdeg(jj) = max(0.1, (pwgi(jj,2)+pwg(jj,2))/pwsat(jj))
338 zkersten(jj) = log10(zsatdeg(jj)) + 1.0
344 zkersten(jj) = (1.0-zfrozen2(jj))*zkersten(jj) + &
345 zfrozen2(jj) *zsatdeg(jj)
349 zcond(jj) = zkersten(jj)*(zcondsat(jj)-pconddryz(jj)) + pconddryz(jj)
353 zhcap(jj) = (1.0-pwsat(jj)) * phcapsoilz(jj) + &
354 pwg(jj,2) * xcl * xrholw + &
355 pwgi(jj,2) * xci * xrholi
359 pcg(jj) = 2.*sqrt(xpi/zcond(jj)/zhcap(jj)/xday)
367 pcg(:) = min( pcg(:), pcgmax )
375 zcf(:) = 2.0 * sqrt( xpi/(xcondwtr*xrholw*xcl*xday) )
378 IF(hsnow_isba ==
'D95' .OR. (hsnow_isba ==
'EBA' .AND. oglacier) )
THEN
381 zlams(:) = xcondi * (psnowrhom(:)/xrholw)**1.885
383 pcs(:) = 2.0 * sqrt( xpi/(zlams(:)*psnowrhom(:)*xci*xday) )
394 pct(:) = 1. / ( (1.-pveg(:))*(1.-ppsng(:)-pffg(:)) / pcg(:) &
395 + pveg(:) *(1.-ppsnv(:)-pffv(:)) / pcv(:) &
404 zff(jj) = pveg(jj)*pffv_nosnow(jj) + (1.-pveg(jj))*pffg_nosnow(jj)
409 pct(jj) = 1. / ( (1.-pveg(jj))*(1.-pffg_nosnow(jj)) / pcg(jj) &
410 + pveg(jj) *(1.-pffv_nosnow(jj)) / pcv(jj) &
411 + zff(jj) / zcf(jj) )
424 zc1sat(:) = pc1sat(:)*sqrt(zwsat1(:)/pwsat(:))
431 WHERE (pwg(:,1) > zwwilt(:))
436 pc1(:) = zc1sat(:) * ( zwsat1(:)/pwg(:,1) )**( 0.5*pbcoef(:) + 1 )
452 IF(hc1dry==
'GB93')
THEN
456 IF (pwg(jj,1) <= zwwilt(jj))
THEN
464 zcw1max(jj) = ( 1.19*zwwilt(jj)-5.09 )*ptg(jj) + (-146.4*zwwilt(jj)+1786.)
468 za(jj) = (-1.815e-2*ptg(jj)+6.41)*zwwilt(jj) &
469 + (6.5e-3*ptg(jj)-1.4)
470 zb(jj) = za(jj)*zwwilt(jj)
471 zdelta(jj) = ( zb(jj)*zb(jj) ) / &
472 ( 2.*log( zcw1max(jj) ) )
474 pc1(jj) = zcw1max(jj)*(1. - 2.*pveg(jj)*( 1.-pveg(jj) )) &
475 *exp( -(pwg(jj,1)-zb(jj))*(pwg(jj,1)-zb(jj)) / &
486 IF (pwg(jj,1) <= zwwilt(jj))
THEN
490 zcw1max(jj) = ( 1.19*zwwilt(jj)-5.09 )*ptg(jj) + (-146.4*zwwilt(jj)+1786.)
499 zy2(jj) = zc1sat(jj)*(zwsat1(jj)/zwwilt(jj))**( 0.5*pbcoef(jj) + 1)
503 zcw1max(jj) = max(max(zcw1max(jj),zy2(jj)),zy1(jj))
507 zlymy1(jj) = log( zcw1max(jj)/zy1(jj))
508 zza(jj) = - log( zy2(jj)/zy1(jj))
509 zzb(jj) = 2. * zx2(jj) * zlymy1(jj)
510 zdelta(jj) = 4. * (zlymy1(jj)+zza(jj)) * zlymy1(jj) * zx2(jj)**2
512 za(jj) = (-zzb(jj)+sqrt(zdelta(jj))) / (2.*zza(jj))
514 zb(jj) = za(jj)**2 / zlymy1(jj)
517 pc1(jj) = zcw1max(jj) * exp( - (pwg(jj,1)-za(jj))**2 / zb(jj) )
530 IF(hksat==
'SGH' .OR. hksat==
'EXP')
THEN
535 zwg2(jj)=pwg(jj,2)*(pcondsat(jj,2)/pcondsat(jj,1))**(1./(2.*pbcoef(jj)+3))
537 zwg2(:)=max(zwg2(:),xwgmin)
549 pc2(jj) = (pc2ref(jj)*zwg2(jj) / ( zwsat(jj)-zwg2(jj) + 0.01 )) &
550 *(1.0-(pwgi(jj,1)/(pwsat(jj)-xwgmin)))
557 zx(jj) = zwg2(jj)/zwsat(jj)
559 pwgeq(jj) = zwg2(jj) - zwsat(jj)*pacoef(jj) &
560 * zx(jj)** ppcoef(jj) &
561 *( 1.-exp(ppcoef(jj)*8.*log(zx(jj))))
572 IF (lhook) CALL dr_hook(
'SOIL',1,zhook_handle)
subroutine soil(HC1DRY, HSCOND, HSNOW_ISBA, OGLACIER, PSNOWRHOM, PVEG, PCGSAT, PCGMAX, PC1SAT, PC2REF, PACOEF, PPCOEF, PCV, PPSN, PPSNG, PPSNV, PFFG, PFFV, PFF, PCG, PC1, PC2, PWGEQ, PCT, PCS, PFROZEN1, PTG, PWG, PWGI, PHCAPSOILZ, PCONDDRYZ, PCONDSLDZ, PBCOEF, PWSAT, PWWILT, HKSAT, PCONDSAT, PFFG_NOSNOW, PFFV_NOSNOW)