c file name: hydro.for subroutine hydro ( layrsn, bmrslp, bcftcv, bcrlai, bdmres, * bczrtd, bhfwsf, & bszlyd, bsdblk, bhrwc, bhrwcs, bhrwcf, * bhrwcw, bh0cb, & bheaep, bsfsan, bsfsil, bsfcla, bh0cng, * bh0cnp, bhzper, & bhzirr, bhzrun, bwudav, bhrsk, bhtsmx, * bhtsmn, bhrwc0, & daysim, bsfald, bsfalw, bszlyt, bwzdpt, * bwtdmx, bwtdmn, 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 DATE: 09/30/93 c MODIFIED: 12/10/93 c MODIFIED: 07/29/95 * This change was done the initialization routine (HINIT.FORt). * The 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 + + + KEY WORDS + + + c hydrology c + + + ARGUMENT DECLARATIONS + + + integer layrsn real bsfald real bsfalw real bsdblk(*) real bhrwc0(*) real bcftcv real bh0cb(*) real bsfcla(*) real bhzirr real bszlyd(*) real bhzper real bhrwc(*) real bhrwcf(*) real bhrwcs(*) real bhrwcw(*) real bcrlai real bheaep(*) real bdmres real bczrtd real rn real bhzrun real bsfsan(*) real bhrsk(*) real bsfsil(*) real bh0cng, bh0cnp real bmrslp c real snmlt real temp0 real bwudav real bhfwsf real bszlyt(*) real bwzdpt real bwtdmx, bwtdmn real bhzwid c + + + ARGUMENT DEFINITIONS + + + c amzele - Average elevation of site (m) c bsdblk - Soil bulk density (Mg/m^3) c bhrwc0 - Water content of the soil surface - darcy supplies c bcftcv - Daily percent crop canopy (decimal) c bh0cb - Power of campbell's water release curve model (?) c bhzirr - Daily irrigation (mm) c bszlyd - Depth to bottom of soil layers (passed as mm converted c to m in hinit) c bhzper - Daily deep percolation (mm/day) c layrsn - Number of soil layers used in simulation c bhrwc - Soil water content (mg/mg) c bhrwcf - Soil water content at fc and wp (mg/mg) c bhrwcs - Soil water content at saturation (mg/mg) c bhrwcw - Soil water content at wp (mg/mg) c bcrlai - Plant leaf area index (unitless) c bheaep - Soil air entry potential (j/kg) c bdmres - Plant residues on the soil surface (kg/ha) c bczrtd - Plant root depth (m) c bhzrun - Daily surface runoff (mm/day) c bsfsan - Mass fraction sand, silt and clay c bszlyt - Soil layer thickness (mm) c bhrsk - Saturated soil hydraulic conductivity (m/s) c bh0cng - Scs curve number for class ii antecedent soil c bh0cnp - moisture conditions under good and poor c hydrologic conditions c bmrslp - Average slope of subregion (mm/mm) c snmlt - Snow melt c temp0 - Local temporary variable used to store water c content of the snow (mm) c bwudav - Wind speed (m/s) c bhfwsf - Plant growth water stress factor (unitless) c + + + COMMON BLOCKS + + + *$noereference include 'p1werm.inc' include 'p1unconv.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1flag.inc' include 'w1clig.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' *$reference c + + + LOCAL VARIABLES + + + integer day, mo, yr integer daysim, l integer idoy integer theta_check c integer subr c integer sunny real rise real bhtsmx(mnsz) real bhtsmn(mnsz) real ts1avg real snwci real lswc save lswc, snwci data lswc /0.0/ data snwci /0.0/ c + + + LOCAL DEFINITIONS + + + c idoy - Day of year c daysim - Current day number of simulation c subr - Current subregion c sunny - Function to determine time of sunrise c l - Loop index of soil layers c bhtsmx - Maximum soil temperature for each layer c bhtsmn - Minimum soil temperature for each layer c ts1avg - Average soil temperature for the surface layer c snwci - Initial snow water content (mm) c + + + SUBROUTINES CALLED + + + c hinit c heat c snomlt c store c et c darcy c transp c caldat C calrwc c + + + FUNCTION DECLARATIONS + + + integer dayear real dawn real radnet real calrwc c + + + DATA INITIALIZATIONS + + + if (daysim .eq. 1) snwci = 0.0 c Calculate hour of sunrise call caldat(-1, day, mo, yr) idoy = dayear(day, mo, yr) rise = dawn(amalat, amalon, idoy) c Initialize others C23456789*123456789*123456789*123456789*123456789*123456789*123456789*1234 call hinit (layrsn, bsdblk, bhrwc, bhrwcs, bhrwcf, * bhrwcw, bh0cb, bheaep, bhrsk, bsfsan, bsfsil, bsfcla, * 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('-')) 2080 format (1x,82('-')/ * 1x,'Sub- Date',t17,'----',t24,'Pot',t31,'----',t38, & '--- Act ---',t52,'Deep', * t80,'Water',t87,'Water',t94,'Check',t101,'Root',t108,'Root', * t115,'Snow'/, * 1x,'Region',t17,'Total', t24,'Evap', t31, 'Trans',t38, 'Evap', * t45,'Trans', t52, 'Perc', t59, 'Irrig', t66, 'Precip', t73, * 'Runoff',t80, 'Soil',t87,'Stress',t101,'Depth',t108,'Water' * ,t115,'Water'/,t17,'*hzetp',t24,'*hzep',t31,'*hzptp',t38 & ,'*hzea',t45,'*hzpta',t52,'*hzper',t59,'*hzirr',t66,'dh2o', * t73,'*hzrun' & ,t80,'swc',t87,'*hfwsf',t101,'*czrtd',t115,'*snwci' & ,t122,'theta'/, * /, t16,34('-'),'mm',33('-')/1x, 100('-')) 2090 format(1x,3(i2,1x),i4,1x,15f7.2,10f7.3) c + + + END SPECIFICATIONS + + + c Open daily output data file (hourly output file is opened in darcy) if ((am0ifl .eqv. .true.) .and.((am0hfl .eq. 1).or.(am0hfl .eq. 3) & .or. (am0hfl .eq. 5) .or. (am0hfl .eq. 7))) then open(14,file= 'hydro.out', & access='sequential', status='unknown') c Echo print of input soil data write(14,2009) write(14,2010) am0csr write(14,2020) do 130 l=1,layrsn write(14,2030) l,bszlyd(l),theta(l),thetas(l),thetaf(l), & thetaw(l),bh0cb(l),bheaep(l),bhrsk(l),bsdblk(l) 130 continue write(14,2040) write(14,2050) theta(0) write(14,2060) swci write(14,2070) write(14,2080) end if c Convert global to net radiation ts1avg = (bhtsmx(1) + bhtsmn(1)) / 2. rn = radnet( bcrlai, aweirr, snwc, awtdmx, awtdmn, amalat, bsfalw, & bsfald, idoy, awtdpt, bwzdpt ) 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 = bwzdpt if ( ((bwtdmn+bwtdmx)/2) .le. 0.0 ) then snwc = snwc + dh2o dh2o = 0.0 endif if ( snwc .gt. 0.0 .and. awtdmx.gt.0) then snmlt = 4.57 * awtdmx if (snmlt .gt. snwc) snmlt = snwc snwc = snwc - snmlt dh2o = dh2o + snmlt endif snwci = snwc c Calculate the amount of water infiltrated into the soil c profile by subtracting surface runoff from the cumulative c amount of rainfall and snow melt. c If measured daily surface runoff is not readily available, c then daily surface runoff is estimated by the scs curve c number method using subroutine scsq. if ( dh2o .lt. 0.0 ) then dh2o= 0.0 endif if ( bh0cnp .eq. 0.0 ) then bhzrun = 0.0 else bhzrun = scsq(dh2o,bh0cnp,bh0cng,bcftcv,bmrslp, & theta(1),thetaf(1)) endif dinf = dh2o - bhzrun c add effective irrigation water dinf = dinf + bhzirr c check that no soil water contents are greater than saturated c this is included until the adding of water and bulk density c reduction can be done simultaneously theta_check = 0 do l=1,layrsn if( theta(l).gt.thetes(l) ) theta_check = 1 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. bhzwid = 0.0 if( (dinf.gt.0.0).or.(theta_check.eq.1) ) then if( theta_check.eq.1 ) write(*,*) 'theta_check' call store(layrsn, bheaep, bh0cb, bszlyt, bhzwid) endif c Calculate potential evapotranspiration using subroutine et call et(rn, bwudav, amzele, bwtdmx, bwtdmn) if ( ahzetp .le. 0.0 ) then ahzep = 0.0 ahzptp = 0.0 ahzetp = 0.0 else c Calculate potential bare soil evaporation ahzep = ahzetp * exp( -0.398 * bcrlai ) !h-35 c Calculate potential plant transpiration if( bcrlai.le.0.09 ) then ahzptp = 0.0 else ahzptp = ahzetp*(-0.21+0.7*((min(bcrlai,3.0))**0.5)) endif 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. temp0 = snwc if ( snwc .gt. 0.0 ) then snwc = snwc - ahzep / 2 ! /2 added per ffox if ( snwc .lt. 0.0 ) then ahzep = -snwc * 2 snwc = 0.0 else ahzep = 0.0 end if 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) call darcy(layrsn, bszlyd, bheaep, bhzper, rise, bhrsk, bh0cb, * rn, bhrwc0, bszlyt, daysim, bsdblk(1)) call timer(TIMDARC,TIMSTOP) call timer(TIMHYDR,TIMSTART) c following darcy, check total et against reduced soil surface ET 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 call heat( layrsn, bszlyd, theta, bsfcla, bsdblk, & bwtdmn, bwtdmx, bhtsmx, bhtsmn, bszlyt) c Print the daily soil water balance results to the output file: c hydro.out swc = dot_product(theta(1:layrsn),bszlyt(1:layrsn)) if ((am0hfl .eq. 1) .or. (am0hfl .eq. 3) .or. (am0hfl .eq. 5) .or. & (am0hfl .eq.7)) then call caldat(-1,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) write(14,2090) am0csr, day, mo, yr, ahzetp, ahzep, ahzptp, & ahzea, ahzpta, bhzper, bhzirr, dh2o, bhzrun, swc, bhfwsf, * lswc+dh2o+bhzirr-swc-ahzea-ahzpta-bhzper-bhzrun, * bczrtd, calrwc(bczrtd*mtomm,bszlyt,theta,thetaw,layrsn), * snwci,theta(1:layrsn) lswc = swc end if c c close ( unit = 14 ) 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