c$Author: fredfox $ c$Date: 2001-11-21 22:38:05 $ c$Revision: 1.25 $ c$Source: /weru/cvs/weps/weps.src/hydro/hinit.for,v $ subroutine hinit(layrsn, bsdblk, bsdblk0, bsdpart, bsdwblk, & bhrwc, bhrwcs, bhrwcf, bhrwcw, bhrwca, & bh0cb, bheaep, bhrsk, & bsfsan, bsfsil, bsfcla, bsfom, bszlyd, bszlyt) c + + + PURPOSE + + + c This subroutine controls the initialization of the HYDROLOGY c sumbodel of WEPS. The program initializes the depth variables c of the soil simulation layers and converts the soil water c content variables from mass basis to volume basis. c DATE: 09/16/93 c MODIFIED: 12/13/93 c MODIFIED: 07/28/95 c MODIFIED: 07/29/95 * This change was done to determine the average of soil * properties from the first simulation layer (10 mm) and the second * uppermost simulation layer (40 mm). The average for the new * uppermost simulation layer for the HYDROLOGY submodel (50 mm thick) * Using 50 mm as the thickness of the uppermost simulation layer for the * HYDROLOGY submodel will increase the speed of simulation and reduce the * the potential for errors. * * c + + + KEYWORDS + + + c initialization, hydrology c + + + ARGUMENT DECLARATIONS + + + integer layrsn real bsdblk(*) real bsdblk0(*) real bsdpart(*) real bsdwblk(*) real bhrwc(*) real bhrwcs(*) real bhrwcf(*) real bhrwcw(*) real bhrwca(*) real bh0cb(*) real bheaep(*) real bhrsk(*) real bsfsan(*) real bsfsil(*) real bsfcla(*) real bsfom(*) real bszlyd(*) real bszlyt(*) c + + + ARGUMENT DEFINITIONS + + + c layrsn - Number of soil layers used in simulation c bsdblk - Soil bulk density (Mg/m^3) c bsdblk0 - Previous day soil bulk density (Mg/m^3) c bsdpart - Soil particle density (Mg/m^3) c bsdwblk - Soil wet bulk density (Mg/m^3) c bhrwc - Soil water content (mg/mg) c bhrwcs - Soil water content at saturation (mg/mg) c bhrwcf - Soil water content at field capacity (mg/mg) c bhrwcw - Soil water content at wilting point (mg/mg) c bh0cb - Power of campbell's water release curve model (unitless) c bheaep - Soil air entry potential (j/kg) c bhrsk - Saturated hydraulic conductivity (m/s) c bsfsan - Sand fractions c bsfsil - Silt fractions c bsfcla - Clay fractions c bszlyd - depth to the bottom of soil layer (mm) c bszlyt - soil layer thickness (mm) c + + + COMMON BLOCKS + + + include 'p1werm.inc' include 'p1unconv.inc' c + + + LOCAL COMMON BLOCKS + + + include 'hydro/htheta.inc' include 'hydro/hpsd.inc' include 'hydro/clayomprop.inc' c + + + LOCAL VARIABLES + + + integer k real potes(mnsz), temp, temp1, temp2, temp3 c + + + LOCAL DEFINITION + + + c potes - Air entry potential at a std. bsdblk of 1.3 Mg/m^3 c + + + SUBROUTINES CALLED + + + c psd c extra c + + + FUNCTION DECLARATIONS + + + real waterk real calctht0 real volwatadsorb real volwat_matpot_bc c + + + END SPECIFICATIONS + + + c initialize various soil layer references c adjust hydro parameters for a change in bulk density call param_blkden_adj( layrsn, bsdblk, bsdblk0, & bsdpart, bhrwcf, bhrwcw, bhrwca, & bsfcla, bsfom, & bh0cb, bheaep, bhrsk ) c convert soil water contents from mass basis to volume basis do k=1,layrsn theta(k) = bhrwc(k) * bsdblk(k) thetas(k) = 1 - bsdblk(k) / bsdpart(k) ! saturation thetes(k) = thetas(k) * 0.883 ! reduced saturation content thetaf(k) = bhrwcf(k) * bsdblk(k) ! field capacity thetaw(k) = bhrwcw(k) * bsdblk(k) ! wilting point c use theta corresponding to 80% relhum in soil for thetar temp = bsdblk(k)*1000.0 !convert Mg/m^3 to kg/m^3 thetar(k) = volwatadsorb( temp, bsfcla(k), bsfom(k), & claygrav80rh, orggrav80rh ) call propsaxt(bsfsan(k), bsfcla(k), & temp, temp1, temp2 ) !thetas, thetaf, thetaw temp3 = (1.0-temp) * bsdpart(k) !bulk density c write(*,1000) k,bszlyt(k), c & bsfsan(k),bsfcla(k),bsfom(k), c & bsdblk(k),thetas(k), c & thetaf(k),thetaw(k), c & temp3, temp, !bulkden, sat vol c & temp1, !field vol c & temp2 !wilt vol c used with output for soil file screening c 1000 format(i3,f7.0,20f7.4) c write(*,*) 'hinit:',k,bh0cb(k),bheaep(k),bhrsk(k),thetas(k), c & thetaf(k),thetaw(k),thetar(k) c c Campbell functions c call psd( bsfsan(k), bsfsil(k), bsfcla(k), gmd(k), gsd(k) ) c potes(k) = -0.2 * gmd(k)**(-0.5) !H-77 c bh0cb(k) = -2. * potes(k) + 0.2 * gsd(k) !H-78 c bheaep(k) = potes(k)*(bsdblk(k)/1.3)**(0.67*bh0cb(k)) !H-79 c c reverse calculation of field capacity and wilting point c thetar(k) = 0.0 c temp = -33.33 c temp1 = 1.0 / bh0cb(k) c thetaf(k) = volwat_matpot_bc(temp, thetar(k), thetas(k), c & bheaep(k), temp1) c temp = -1500.0 c thetaw(k) = volwat_matpot_bc(temp, thetar(k), thetas(k), c & bheaep(k), temp1) c c write(*,*) 'hinit:',k,bsfsan(k),bsfsil(k),bsfcla(k),bsfom(k), c & bh0cb(k),bheaep(k),thetas(k), c & thetaf(k),thetaw(k),thetar(k),bhrsk(k) c c this is used with campbell functions as well c bhrsk(k) = waterk(bsdblk(k), bh0cb(k), bsfcla(k), bsfsil(k)) c c write(*,*) 'hinit:',k,bsfsan(k),bsfsil(k),bsfcla(k),bsfom(k), c & bh0cb(k),bheaep(k),thetas(k), c & thetaf(k),thetaw(k),thetar(k),bhrsk(k) end do c swci = sum(wc(1:layrsn)) swci = dot_product(theta(1:layrsn),bszlyt(1:layrsn)) c theta(0) = calctht0(bszlyd, thetes, thetar, theta, c * thetaw, thetaf(1) - thetaw(1), 0.0_8, 0.0_8) !H-64,65,66 theta(0) = theta(1) c write(*,*) 'calctht0: hinit - theta(0) after',theta(0) c call subroutine psd and calculate soil hydraulic parameters do 110 k=1,layrsn if ( bh0cb(k) .eq. -99.9 ) then call psd( bsfsan(k), bsfsil(k), bsfcla(k), gmd(k), * gsd(k) ) potes(k) = -0.2 * gmd(k)**(-0.5) !H-77 bh0cb(k) = -2. * potes(k) + 0.2 * gsd(k) !H-78 if ( bheaep(k) .eq. -99.90 ) then bheaep(k) = potes(k)*(bsdblk(k)/1.3)**(0.67*bh0cb(k)) !H-79 end if end if if ( bhrsk(k) .eq. -99.90 ) then bhrsk(k) = waterk(bsdblk(k), bh0cb(k), bsfcla(k), * bsfsil(k)) end if C *** write (*,*) ' theta(k), etc. ', k, theta(k),thetas(k),-bh0cb(k) 110 continue return end * *