c$Author: fredfox $ c$Date: 2001-08-15 18:28:19 $ c$Revision: 1.55.2.3 $ c$Source: /weru/cvs/weps/weps.src/hydro/hydro.for,v $ c file name: hydro.for subroutine hydro ( layrsn, bmrslp, bcftcv, bcrlai, bdmres, * bczrtd, bhfwsf, & bszlyd, bsdblk, bsdpart, bsdwblk, & bhrwc, bhrwcs, bhrwcf, * bhrwcw, bh0cb, bheaep, & bsfsan, bsfsil, bsfcla, & bsfom, bhtsav, bbdstm, bbffcv, & bszrro, bszrr, bmzele, & bh0cng, bh0cnp, bhzper, & bhzirr, bhzrun, bhrsk, & bhtsmx, bhtsmn, bhrwc0, & daysim, bsfald, bsfalw, bszlyt, & bwzdpt, bwdurpt, bwpeaktpt, * bwtdmxprev, bwtdmn, bwtdmx, bwtdmnnext, & bwtdpt, bweirr, bwudav, bhzwid ) c + + + PURPOSE + + + c This subroutine is the main (supervisory) program for the c HYDROLOGY submodel. The subroutine controls the calling of the c major subprograms of the HYDROLOGY submodel. c + + + KEY WORDS + + + c hydrology c + + + ARGUMENT DECLARATIONS + + + integer layrsn real bmrslp, bcftcv, bcrlai, bdmres, & bczrtd, bhfwsf, & bszlyd(*), bsdblk(*), bsdpart(*), bsdwblk(*), & bhrwc(*), bhrwcs(*), bhrwcf(*), & bhrwcw(*), bh0cb(*), bheaep(*), & bsfsan(*), bsfsil(*), bsfcla(*), & bsfom(*), bhtsav(*), bbdstm, bbffcv, & bszrro, bszrr, bmzele, & bh0cng, bh0cnp, bhzper, & bhzirr, bhzrun, bhrsk(*), & bhtsmx(*), bhtsmn(*), bhrwc0(*) integer daysim real bsfald, bsfalw, bszlyt(*) real bwzdpt, bwdurpt, bwpeaktpt real bwtdmxprev, bwtdmx, bwtdmn, bwtdmnnext real bwtdpt, bweirr, bwudav, bhzwid c + + + ARGUMENT DEFINITIONS + + + c layrsn - Number of soil layers used in simulation c bmrslp - Average slope of subregion (mm/mm) c bcftcv - Daily percent crop canopy (decimal) c bcrlai - Plant leaf area index (unitless) c bdmres - Plant residues on the soil surface (kg/ha) c bczrtd - Plant root depth (m) c bhfwsf - Plant growth water stress factor (unitless) c bszlyd - Depth to bottom of soil layers (mm) 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 fc and wp (mg/mg) c bhrwcw - Soil water content at wp (mg/mg) c bh0cb - Power of campbell's water release curve model (?) c bheaep - Soil air entry potential (j/kg) c bsfsan - Fraction of soil mineral which is sand c bsfsil - Fraction of soil mineral which is silt c bsfcla - Fraction of soil mineral which is clay c bsfom - Fraction of total soil mass which is organic c bhtsav - daily average soil temperature (C) c bbdstm - total number of stems (#/m^2) c bszrro - original random roughness height, after tillage, mm c bszrr - Allmaras random roughness parameter (mm) c bmzele - Average site elevation (m) c bh0cng - Scs curve number for class ii antecedent soil c bh0cnp - moisture conditions under good and poor c hydrologic conditions c bhzper - Daily deep percolation (mm/day) c bhzirr - Daily irrigation (mm) c bhzrun - Daily surface runoff (mm/day) c bhrsk - Saturated soil hydraulic conductivity (m/s) c bhtsmx - Maximum soil temperature for each layer c bhtsmn - Minimum soil temperature for each layer c bhrwc0 - Hourly water content of the soil surface - darcy supplies c daysim - day of the simulation c bsfald - dry soil albedo c bsfalw - wet soil albedo c bszlyt - Soil layer thickness (mm) c bwzdpt - Daily precipitation (mm) c bwdurpt - Duration of Daily precipitation (hours) c bwpeaktpt - Normalized time to peak of Daily precipitation (time to peak/duration) c bwtdmxprev - previous day maximum daily air temperature (deg C) c bwtdmn - Minimum daily air temperature (deg C) c bwtdmx - Maximum daily air temperature (deg C) c bwtdmnnext - Following day minimum daily air temperature (deg C) c bwtdpt - Daily dewpoint temperature (C) c bweirr - Daily global total solar radiation (MJ/m^2) c bwudav - Daily average wind speed at 10 meters (m/s) c bhzwid - Water infiltration depth (mm) c + + + COMMON BLOCKS + + + include 'p1werm.inc' include 'p1unconv.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1flag.inc' include 'h1et.inc' include 'h1hydro.inc' include 'h1db1.inc' include 'timer.fi' c + + + LOCAL COMMON BLOCKS + + + include 'hydro/hfdecl.inc' include 'hydro/htheta.inc' include 'hydro/hh2o.inc' c + + + LOCAL VARIABLES + + + integer day, mo, yr integer l integer idoy integer theta_check real rise, daylength real ts1avg real availwat(layrsn) real rn integer numeq c + + + LOCAL DEFINITIONS + + + c idoy - Day of year c daysim - Current day number of simulation c l - Loop index of soil layers c rise - time of sunrise adjusted to local time (hours) c daylength - length of day (hours) c ts1avg - Average soil temperature for the surface layer c rn - net radiation from function radnet c numeq - number of equations to be solved by darcy routines c + + + SUBROUTINES CALLED + + + c hinit c heat c snomlt c store c et c darcy c transp c caldat c + + + FUNCTION DECLARATIONS + + + integer dayear real dawn real daylen real radnet real calrwc real calrwcap real availwc real plant_wat_t c + + + DATA INITIALIZATIONS + + + c Calculate hour of sunrise call caldatw(day, mo, yr) idoy = dayear(day, mo, yr) rise = dawn(amalat, amalon, idoy) daylength = daylen(amalat, idoy) c write(*,*)'hydro:',daysim,rise,daylength c write(*,*)'hydro:rain:',daysim,bwzdpt,bwdurpt,bwpeaktpt c Initialize others if(daysim.eq.1) lswc = 0.0 if(daysim.eq.1) snwci = 0.0 call hinit (layrsn, bsdblk, bsdpart, bsdwblk, bhrwc, bhrwcs, * bhrwcf, bhrwcw, bh0cb, bheaep, bhrsk, bsfsan, bsfsil, bsfcla, * bsfom, bszlyd, bszlyt) c + + + OUTPUT FORMATS + + + 2000 format(' ') 2009 format(' Daily HYDROLOGY output ') 2010 format (/,1x,18x,5('*'),' soil data - Subregion #',i4,3x, & 5('*'),16x) 2020 format (1x,79('-')/' soil depth',t16,'initial',t25,'saturated', & t36,'field',t45,'wilting',t54,'bh0cb',t60,'air',t69,'sat.',t76, & 'bulk'/' layer',t17,'water',t27,'water',t35, & 'capacity',t46,'point term entry',t70,'k density'/' no', & t16,'content',t26,'content',t57,'potential'/t9, & '(mm)',3x,15('-'),'m^3/m^3',14('-'),t57,'joules/kg', & t69,'m/s',3x,'mg/m^3',/1x,79('-')) 2030 format (1x,i3,f8.3,f9.3,f10.3,2f9.3,2f7.2,e10.3,f6.2) 2040 format (1x,79('-')) 2050 format(1x,'initial soil wetness at the soil-air interface =', & f6.3,' m^3/m^3') 2060 format(1x,'initial total amount of soil water in the soil' & ,' profile =',f7.2,' mm'/) 2070 format (/1x,19x,5('*'),' daily soil water balance data ', & 5('*')) C *** 2080 format (1x,82('-')/1x,'sr date',t17,'etp',t23,'ep',t29,'tp',t36 C *** & ,'ea',t43,'ta',t47,'bhzper',t54,'bhzirr',t61,'dh2o',t67, C *** * 'bhzrun' C *** & ,t74,'swc',t79,'bhfwsf',' check '/,t16,31('-'),'mm',30('-')/1x, C *** * 82('-')) c 2080 format (1x,82('-')/ c * 1x,'Sub- Date',t17,'----',t24,'Pot',t31,'----',t38, c & '--- Act ---',t52,'Deep', c * t80,'Water',t87,'Water',t94,'Check',t101,'Root',t108,'Root', c * t115,'Snow'/, c * 1x,'Region',t17,'Total', t24,'Evap', t31, 'Trans',t38, 'Evap', c * t45,'Trans', t52, 'Perc', t59, 'Irrig', t66, 'Precip', t73, c * 'Runoff',t80, 'Soil',t87,'Stress',t101,'Depth',t108,'Water' c * ,t115,'Water'/,t17,'*hzetp',t24,'*hzep',t31,'*hzptp',t38 c & ,'*hzea',t45,'*hzpta',t52,'*hzper',t59,'*hzirr',t66,'dh2o', c * t73,'*hzrun' c & ,t80,'swc',t87,'*hfwsf',t101,'*czrtd',t115,'*snwci' c & ,t122,'theta'/, c * /, t16,34('-'),'mm',33('-')/1x, 100('-')) 2080 format(' am0csr day mo yr ahzetp ahzep ahzptp ahzea ahzpta bhzper & bhzirr dh2o bhzrun swc bhfwsf check bczrtd plntwc plntcap snwci & theta') 2090 format(1x,i5,1x,2(i2,1x),i4,1x,16f7.2,400f7.3) c + + + END SPECIFICATIONS + + + c write headers and inital values to hydro.out if ((am0ifl .eqv. .true.) .and.((am0hfl .eq. 1).or.(am0hfl .eq. 3) & .or. (am0hfl .eq. 5) .or. (am0hfl .eq. 7))) then c Echo print of input soil data c write(14,2009) c write(14,2010) am0csr c write(14,2020) c do 130 l=1,layrsn c write(14,2030) l,bszlyd(l),theta(l),thetas(l),thetaf(l), c & thetaw(l),bh0cb(l),bheaep(l),bhrsk(l),bsdblk(l) c 130 continue c write(14,2040) c write(14,2050) theta(0) c write(14,2060) swci c write(14,2070) write(14,2080) end if c write(*,*) 'hydro:total 500mm', c & plant_wat_t( 0.0, 500.0, thetaf, thetaw, bszlyt, layrsn ), c & plant_wat_t( 500.0, 1000.0, thetaf, thetaw, bszlyt, layrsn ), c & plant_wat_t( 1000.0, 1500.0, thetaf, thetaw, bszlyt, layrsn ) c Determine soil temperature for each layer (before darcy) call heat( layrsn, bszlyd, theta, bsfcla, bsdblk, & bwtdmn, bwtdmx, bhtsmx, bhtsmn, bszlyt) c write(*,*)'hydro:t',daysim,bwtdmn,bwtdmx,bhtsmn(2),bhtsmx(2) c Partition daily precipitation between rainfall and snowfall. c estimate daily snow melt, and add the melted snow if any to c the daily precipitation using subroutine snomlt snwc = snwci dh2o = 0.0 if( bwzdpt.gt.0.0 ) then if ( ((bwtdmn+bwtdmx)/2) .le. 0.0 ) then snwc = snwc + bwzdpt else dh2o = dh2o + bwzdpt endif endif c add effective irrigation water to rainfall if( bhzirr.gt.0.0 ) then if ( ((bwtdmn+bwtdmx)/2) .le. 0.0 ) then snwc = snwc + bhzirr else dh2o = dh2o + bhzirr endif c when irrigation water is added, it must have a duration (hours) c since rainfall may already have a duration, it needs to be added bwdurpt = max(bwdurpt,12.0) !we have no entry for irrigation duration yet bwpeaktpt = 0.5 endif c add snowmelt to surface water for infiltration if ( (snwc.gt.0.0).and.(bwtdmx.gt.0.0) ) then snmlt = 4.57 * bwtdmx if (snmlt .ge. snwc) then snmlt = snwc snwc = 0.0 else snwc = snwc - snmlt end if dh2o = dh2o + snmlt bwdurpt = max(bwdurpt,6.0)!we have no entry for snowmelt duration yet bwpeaktpt = 0.5 endif snwci = snwc c check that no soil water contents are greater than saturation c this is included until the adding of water and bulk density c reduction can be done simultaneously c *** temporarily disabled c *** theta_check = 0 c *** do l=1,layrsn c *** if( theta(l).gt.thetas(l) ) theta_check = 1 c *** end do c Distribute the total amount of water available for c infiltration throughout the simulation layers of the soil c profile using subroutine store. c *** bhzwid = 0.0 c *** if( theta_check.eq.1 ) then c if( theta_check.eq.1 ) write(*,*) 'theta_check' c *** call store(layrsn, bheaep, bh0cb, bszlyt, bhzwid) c *** endif c Convert global to net radiation ts1avg = (bhtsmx(1) + bhtsmn(1)) / 2. rn = radnet( bcrlai, bweirr, snwc, bwtdmx, bwtdmn, amalat, & bsfalw, bsfald, idoy, bwtdpt, bwzdpt ) c Calculate potential evapotranspiration using subroutine et call et(rn, bwudav, bmzele, bwtdmx, bwtdmn) if ( ahzetp .le. 0.0 ) then ahzep = 0.0 ahzptp = 0.0 ahzetp = 0.0 else c Calculate potential plant transpiration ahzptp=ahzetp*max(0.0,min(1.0,(-0.21+0.7*(bcrlai**0.5)))) ahzep = ahzetp - ahzptp c If snow is present, it is evaporated at the potential soil c evaporation rate. the potential soil evaporation is then c adjusted by accounting for snow evaporation. if ( snwc .gt. 0.0 ) then if( ahzep.gt.snwc ) then ahzep = ahzep - snwc snwc = 0.0 else snwc = snwc - ahzep ahzep = 0.0 endif end if c reinitialize snow water content snwci = snwc c If plant residue is present, then reduce the potential soil c evaporation on the basis of the amount of plant residues on c the soil surface. if ( ahzep .gt. 0.0 ) then if ( bdmres .gt. 0.0 ) then ahzep = ahzep * exp(-0.64 * bdmres) !h-37 bdmres in kg/m^2 end if end if endif c update global snow variables ahzsno = snwc ahzsnd(am0csr) = snwc * 5.0 !this is a gross approximation, look for better c the following 1 line added by Retta on 7/19/1996 after noticing the c fact that transpiration was being calculated during the dormancy c period of w. wheat. the problem was pointed out to E. Skidmore and c A. Durar. c set potential plant transpiration to 0 during dormancy if (am0drmfl.eq.1) ahzptp = 0.0 c Calculate soil water redistribution using subroutine darcy call timer(TIMHYDR,TIMSTOP) call timer(TIMDARC,TIMSTART) numeq = layrsn + 6 call darcy( daysim, numeq, bszlyt, bszlyd, theta, bsdblk, & thetas, thetar, bhrsk, bheaep, bh0cb, & bsfcla, bsfom, bhtsav, & bwtdmxprev, bwtdmn, bwtdmx, bwtdmnnext, bwtdpt, & rise, daylength, ahzep, dh2o, bwdurpt, bwpeaktpt, & bbdstm, bbffcv, bszrro, bszrr, bmzele, bhrwc0, & ahzea, bhzper, bhzrun, bhzwid ) call timer(TIMDARC,TIMSTOP) call timer(TIMHYDR,TIMSTART) c following darcy, check total et against reduced soil surface ET c ahzptp = min(ahzptp, ahzetp - ahzea) c Calculate actual plant transpiration using subroutine transp c and determine the water stress factor if ( ahzptp .gt. 0.0 ) then call transp(layrsn, bszlyd, bszlyt, bczrtd, bhfwsf) else ahzpta = 0.0 bhfwsf = 1.0 end if c Calculate actual evapotranspiration ahzeta = ahzea + ahzpta c Convert water from volume back to mass basis do l=1,layrsn bhrwc(l) = theta(l) / bsdblk(l) end do c Determine soil temperature for each layer (again with new water content) call heat( layrsn, bszlyd, theta, bsfcla, bsdblk, & bwtdmn, bwtdmx, bhtsmx, bhtsmn, bszlyt) c Print the daily soil water balance results to hydro.out if ((am0hfl .eq. 1).or.(am0hfl .eq. 3).or.(am0hfl .eq. 5).or. & (am0hfl .eq.7)) then do l=1,layrsn availwat(l) = availwc(theta(l),thetaw(l),thetaf(l)) end do swc = dot_product(theta(1:layrsn),bszlyt(1:layrsn)) call caldatw(day,mo,yr) if (lswc.eq.0.0) lswc = -(dh2o+bhzirr-swc-ahzea * -ahzpta-bhzper-bhzrun) if ((am0csr .eq. 1) .and. (nsubr .gt. 1)) write(14,2000) c 2090 format(1x,i5,1x,2(i2,1x),i4,1x,16f7.2,40f7.3) write(14,2090) daysim, day, mo, yr, ahzetp, ahzep, ahzptp, & ahzea, ahzpta, bhzper, bhzirr, dh2o, bhzrun, swc, bhfwsf, * swc-lswc, c * lswc+dh2o+bhzirr-swc-ahzea-ahzpta-bhzper-bhzrun, * bczrtd,calrwc(bczrtd*mtomm,bszlyt,theta,thetaw,layrsn), * calrwcap(bczrtd*mtomm,bszlyt,thetaf,thetaw,layrsn), * snwci,theta(1:layrsn),thetaf(1:layrsn),thetaw(1:layrsn), * availwat(1:layrsn) lswc = swc end if 110 continue return end C********************************************************************wjr real function calrwc(bczrtd,bszlyt,btheta,bthetaw,layrsn) real bczrtd !root depth (mm) real bszlyt(*) !layer thicknesses (mm) real btheta(0:*) !vol water content (mm/mm) real bthetaw(*) !vol water wilting point (mm/mm) integer layrsn !number of layers used integer daysim real depth real tdepth real tswc integer idx tswc = 0.0 tdepth = bczrtd do idx = 1,layrsn depth = min(tdepth, bszlyt(idx)) tswc = tswc + depth * (btheta(idx) - bthetaw(idx)) tdepth = tdepth - depth end do calrwc = tswc end C********************************************************************wjr real function calrwcap(bczrtd,bszlyt,btheta,bthetaw,layrsn) real bczrtd !root depth (mm) real bszlyt(*) !layer thicknesses (mm) real btheta(*) !vol water content (mm/mm) real bthetaw(*) !vol water wilting point (mm/mm) integer layrsn !number of layers used integer daysim real depth real tdepth real tswc integer idx tswc = 0.0 tdepth = bczrtd do idx = 1,layrsn depth = min(tdepth, bszlyt(idx)) tswc = tswc + depth * (btheta(idx) - bthetaw(idx)) tdepth = tdepth - depth end do calrwcap = tswc end