c$Author: fredfox $ c$Date: 2001-09-27 16:15:48 $ c$Revision: 1.23.2.4 $ c$Source: /weru/cvs/weps/weps.src/hydro/hinit.for,v $ subroutine hinit (layrsn, bsdblk, bsdpart, bsdwblk, bhrwc, & bhrwcs, bhrwcf, bhrwcw, 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 bsdpart(*) real bsdwblk(*) real bhrwc(*) real bhrwcs(*) real bhrwcf(*) real bhrwcw(*) 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 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 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 convert soil water contents from mass basis to volume basis do 105 k=1,layrsn theta(k) = bhrwc(k) * bsdblk(k) thetaf(k) = bhrwcf(k) * bsdblk(k) ! field capacity thetaw(k) = bhrwcw(k) * bsdblk(k) ! wilting point thetas(k) = 1 - bsdblk(k) / bsdpart(k) ! saturation thetes(k) = thetas(k) * 0.883 ! reduced saturation content c write(*,*) 'hinit:',k,bh0cb(k),bheaep(k),bhrsk(k),thetas(k), c & thetaf(k),thetaw(k),thetar(k) 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 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 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 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 ) c error check and adjustments to prevent numerical problems c with curve fit if( thetas(k).le.thetaf(k) ) then write(*,*) 'saturation less than field capacity, layer',k stop else if( thetaf(k).le.thetaw(k) ) then write(*,*) 'field capacity less than wilting point, layer',k stop end if thetar(k) = min( thetar(k), 0.8 * thetaw(k) ) c Calculate air entry and b to match saturated, field c capacity, permanent wilting point and residual water c as calcuated above bh0cb(k) = -(log(33.333)-log(1500.0)) & / (log((thetaf(k)-thetar(k))/(thetas(k)-thetar(k))) & - log((thetaw(k)-thetar(k))/(thetas(k)-thetar(k)))) bheaep(k)=-exp(log(1500.0)+bh0cb(k) & *log((thetaw(k)-thetar(k))/(thetas(k)-thetar(k)))) c this is used with campbell functions as well c bhrsk(k) = waterk(bsdblk(k), bh0cb(k), bsfcla(k), bsfsil(k)) 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) 105 continue 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 * *