program bnmf2v c modify bnmf2u c set up DIRK2 time step, w/ simple jacobian implicit real*8(a-h,o-z) character*7 outfs(1) dimension istr(8) character*10, cp1, cp2, cp3 outfs(1) = 'jn1700f' atr = 1.d0 vnr = 1.d0 vmr = 1.0d0 scr0 = 1.d0 bcr = 1.d0 open(unit=2,file=outfs(1)//'1.dat', status='unknown') r = 10.d0**(0.1d0) c call date_and_time (cp1,cp2,cp3,istr) c stime = 0.001d0*istr(8)+istr(7)+60.d0*(istr(6)+60.d0*(istr(5) c 1 +24.d0*istr(3))) call timer (istime) stime = 0.01d0*dble(istime) c do ibc = 0, 20 do ibc = 0, 0 scr = scr0 c do isc = 0, 20 do isc = 0, 0 call bnm (atr, vnr, vmr, scr, bcr, sofp, soft, outfs(1), + stime, noluck) write (2, '(f7.4, f8.4, g14.6, i3)') scr, bcr, sofp, noluck write (*, '(f7.4, f8.4, g14.6, i3)') scr, bcr, sofp, noluck scr = scr*r end do bcr = bcr*r end do close (2) pause end subroutine bnm (atr, vnr, vmr, scr, bcr, sofp, soft, outfs, strt, + noluck) c Interpolated soil hydraulic parameters and properties c Put rx(i) back in c c c ****************************************************************** c * c * simulation of nitrogen transformation and transport c * c * from land disposal of poultry litter c * c * h. d. scott, d. l. baker & b. a. ibrahim c * department of agronomy c * university of arkansas c * 115 plant science c * fayetteville, ar 72701 c * usa c * c * telephone (501) 575-5740 c * c * fax (501) 575-7465 c * c * version 3.0 c * 1996 c ****************************************************************** implicit real*8(a-h,o-z) c*------ bring in the common blocks to be shared between the subroutines c include 'come4.for' c*********************************************************************** c* file come4.for - MULITYEAR SIMULATION c*********************************************************************** c* c* this common blocks file contains all the variables that need to c* be passed between the different subroutines c* c*********************************************************************** parameter (izm=107, iwm=50) common /pr1/ aa(0:izm), bb(0:izm), cc(0:izm), c(0:izm), ch(0:izm), + ck(0:izm), cv(0:izm), ht(0:izm), p(0:izm), qme(0:izm), r(0:izm), + rdist(0:izm), rkd(0:izm), rx(0:izm), t(0:izm), tdif(0:izm), + th(0:izm), thj(0:izm), thr(0:izm), udl(0:izm), y(0:izm), + v(0:izm), z(0:izm), pm(5), tm(2), wm(2), oft(3), ofp(11), + at(0:izm), bc(0:izm), rdn(0:izm), rex(0:izm), rnt(0:izm), + rob(0:izm), sc(0:izm), ths(0:izm), vm(0:izm), vn(0:izm), + thp(0:izm), dk(0:izm), wv(0:izm) common /pr2/ cl, cn, disp, dt, dz, dzl, dzp, fcu, pl, + pmax, qk, ql, qm, qu, res, sev, wilt, wfin, winit, wk1, wk2, + zl, dzg, tdg, dtt, tl, ztl, prk, a1, a2, qp, dtmin, zr, zrp, + zrj, qup, qe, resp, qn, qnp, qin, qinp, qlp, srxp common /pr4/ lt, nm1, n, np1, inewt, maxit common /ts1/ crain(0:iwm), crdaly(-4:iwm), csno3(0:iwm), + et(0:iwm), fam(0:iwm), fav(0:iwm), pet(0:iwm), + runof(0:iwm), tappn(0:iwm), tempmn(0:iwm), tempmx(0:iwm), + xinfl(0:iwm) common /ts2/ accevp, accpln, ampd, aplr, avpd, blt, cav, cavl, + cev, cev1, csin, csinl, csls, cslr, csom, cson, csonl, dnitr1, + dnitrf, fci, gwn, gwn1, gwn2, gwnh4, gwno3, gww1, gww2, + percen, pln, pln1, plnh4, plnh41, plno3, plno31, plnw, plnw1, + qlnh4, qlno3, rad, sumcn, sumtns, sumyn, tave, tavek, + time, timnw, tinf, tmax, tmaxf, tmin, tminf, tnfin, tninit, + twrite, xkt, myr, mdoy, cslo common /ts3/ kdoy(0:iwm), kyer(0:iwm) common /ts4/ imm, idd, iyr, kday, klm, kyr, lqtr, nday, napp logical qtr, rday, vd common /lo1/ qtr, rday, vd(0:izm) dimension rain(0:iwm), istr(8) logical napday, allok character*7 outfs(1) character*10, cp1, cp2, cp3 integer ifc c character*7 outfs(5) c*------ bring in the initial conditions to be used in the simulation noluck = 0 c call timer(istime) c call date_and_time (cp1,cp2,cp3,istr) c istime = 0.001d0*istr(8)+istr(7)+60.d0*(istr(6)+60.d0*(istr(5) c 1 +24.d0*istr(3))) C C fset.f - file unit definitions C c open(unit=32, file='control.dat', status='unknown') open(unit=20, file='weather.dat', status='unknown') open(unit=15, file='soils.dat', status='unknown') open(unit=16, file='sinit.dat', status='unknown') open(unit=77, file='appinfo.dat', status='unknown') open(unit=1, file='outday.dat', status='unknown') rewind (20) rewind (15) rewind (16) rewind (77) rewind (1) c read (32, '(a7)') outfs(1) c open(unit=12,file=outfs(1)//'1.out', status='unknown') c open(unit=14,file=outfs(1)//'2.out', status='unknown') open(unit=25,file=outfs(1)//'3.out', status='unknown') c open(unit=18,file=outfs(1)//'4.out', status='unknown') c open(unit=7,file=outfs(1)//'5.out', status='unknown') c open(unit=2,file=outfs(1)//'6.out', status='unknown') c open(unit=3,file=outfs(1)//'7.out', status='unknown') c open(unit=4,file=outfs(1)//'8.out', status='unknown') c include 'INITCON.FOR' c**************************************************************************** c* file initcon.f * c**************************************************************************** c* * c* this is the include file for reading input data regarding initial * C* time steps and soil depth steps as well as total soil profile depth * C* and depths of individual soil layers. Also constants related to * C* soil water paramters, bulk density, saturation, and exchangeable * C* coefficients for NH4, nitrification, and denitrification are read * C* for each soil layer. Finally time allowed for moisture infiltration * C* and the daily time cycle are read along with the time between * C* detailed information output. All this data is read from the file * C* soils.dat * c* * c**************************************************************************** c* ----------------- read the input data from a prepared data file read (15,*) dtt, dzz, dzp c write (12,101) dtt, dzz, dzp read (15,*) qm, qk, disp c write (12,105) qm, qk, disp read (15,*) cl c write (12,110) cl read (15,*) zl, pl, wk1 c write (12,165) zl, pl, wk1 read (15, *) dzg, tdg, tl, ztl c write (12,166) dzg, tdg, tl, ztl read (15,*) tinf, twrite c write (12,155) twrite c*------- reset boundaries to coincide with grid cell centers dtt = 6.d0 c dzz = 1.d0 dz = dzz n = (cl+0.50*dz)/dz + 0.5d0 cl = n*dz - 0.5d0*dz lt = (ztl + 0.5d0*dz)/dz + 0.5d0 ztl = lt*dz - 0.5d0*dz nm1 = n-1 np1 = n+1 c*------- set up vertical depth array z(0) = 0.0d0 do i = 1, n zz = i*dz - 0.5d0*dz z(i) = zz end do z(np1) = zl c*------- read in soil hyraulic parameters and properties call idist1 (sc) call idist1 (bc) call idist1 (ths) call idist1 (thr) call idist1 (at) call idist1 (vn) call idist1 (vm) call idist1 (rob) call idist1 (rex) call idist1 (rnt) call idist1 (rdn) c*------- add kludges do i = 0, np1 at(i) = atr*at(i) vn(i) = vnr*vn(i) vm(i) = vmr*vm(i) sc(i) = scr*sc(i) bc(i) = bcr*bc(i) end do c*--------- write out the input parameters in a nicely formatted way 100 format(8f5.3) 101 format(5x,'initial dt, hr =',49x,f10.5/, + 5x,'initial dz, cm =',49x,f10.5/,5x,'printout dz, cm=',49x, + f10.5//) 200 format(2f10.5, i3) 300 format(2i4) 600 format(10e12.4) 500 format(e10.4,6f10.5) 105 format(5x, +'nitrogen uptake rate, microgram-n/cm of root length per hour', +' =',3x,f10.5/, +5x,'michaelis constant =',45x,f10.5/, +5x,'concentration of nh4-n in applied litter, ug/cm**2 hr =',10x, +f10.5,/,5x,'concentration of no3-n in applied litter, ug/cm**2 hr +=',10x,f10.5/,5x,'solute dispersion coefficient, cm**2/hr =',24x, +f10.5/) 110 format(5x,'total length of soil profile, cm =',26x,f10.5/) 165 format (19x, ' Depth of lower boundary (cm), zl =', e15.4/ +,12x, ' Pressure head at zl (cm), pl =', e15.4/ +,12x, ' Darcian mean weight, wk1 =,' e15.4/) 166 format (19x, ' Height of insulating grass (cm), dzg =', e15.4/ +,12x, ' Thermal diffusion of grass (cm**2/h), tdg =', e15.4/ +,12x, ' Isotherm temperature (deg-C), tl =', e15.4/ +,12x, ' Isotherm position (cm), ztl =', e15.4/) 150 format(5x,'duration of waste water application, hrs =',f10.5/, +5x,'schedule of waste water application, I.E. cycle duration +=',f10.5/, +5x,'number of cycles =',i3/) 155 format(5x,'time at which output data is required =',26x,f10.5/) c write (12,160) 160 format(t30,'input data ',///) C INCLUDE 'outhead.for' c* print headers for water and nitrogen balance output files c*------- prepare the mass balance table for soil nitrogen content -- c*------- 14 --> *2.out c*------- prepare the mass balance table for soil water content c*------- 25 --> *3.out c*------- 18 --> *4.out litter compartment changes c write (18, 218) 218 format(40x, 'Mineralization & Volatilization Summary') c write (18,*) c write (18,219) 219 format( +' sod eod minerlzd sod eod leached + leached leachate volatil. since removed daily'/, +' organic organic per inorgan inorgan to + to concen- per last to soil avg'/, +' N N day N N soil + runoff tration day rain OM air'/, +'yr day napp csonl cson ampd csinl csin csls + cslo c(0) avpd cav csom temp'/, +' kg/ha kg/ha kg/ha kg/ha kg/ha kg/ha + kg/ha ug/ml kg/ha kg/ha kg/ha deg-C'//, +'================================================================= +====================================================='//) c*------- 7 --> *5.out header for Dr. Cochran's data c include 'INITPRO.FOR' ***************************************************************************** * * * file INITPRO.F * * * c**************************************************************************** c* * c* this file initializes all the variables that are used * c* in the main program. the necessary information for * c* the poultry litter inputs are also contained here. * c* the initial values for soil water pressure, nh4-n * c* concentration, and no3-n concentration will be read * c* from file "sinit.dat" and a linear interpolation function * c* is used to generate the initial profile distribution of * c* these variables. * c* * c**************************************************************************** c* -------- initialization of variables to be used in the ----------- c* -------------------- different subroutines ----------------------- c zr = 0.d0 c zrp = 0.d0 c zrj = 0.d0 dtmin = 1.d-1 maxit = 48 prk = 0.4d0 a1 = 0.5d0/(1.d0-prk) a2 = (1.d0-2.d0*prk)*a1 c prk = 0.d0 c a1 = 0.d0 c a2 = 1.d0 wk1 = 0.5d0 c zl = 1000.d0 c ztl = 171.d0 twrite = 24.d0 tinf = 10.d0 xtinf = tinf time = 0.0d0 sofp = 0.d0 soft = 0.d0 plno3 = 0.0d0 plnh4 = 0.0d0 plno31 = 0.0d0 plnh41 = 0.0d0 dnitrf = 0.0d0 plnw = 0.0d0 plnw1 = 0.0d0 gww1 = 0.0d0 gwn1 = 0.0d0 gww2 = 0.0d0 gwn2 = 0.0d0 gwnh4 = 0.0d0 c gwnh41 = 0.d0 gwno3 = 0.0d0 c gwno31 = 0.d0 pln1 = 0.0d0 cev = 0.0d0 cev1 = 0.0d0 dnitr1 = dnitrf dt = dtt dz = dzz dzl = zl - cl c timel = 0.0d0 timnw = twrite wk2 = 1.0d0 - wk1 pmax = 7.62d0 wilt = -15000.d0 klm = 1 kyer(0) = 0 kdoy(0) = 0 do i = 1, 3 oft(i) = 0.d0 end do do i = 1, 9 ofp(i) = 0.d0 end do do i = 1, iwm crain(i) = 0.d0 fam(i) = 0.d0 fav(i) = 0.d0 tappn(i) = 0.d0 end do cav = 0.d0 cavl = 0.d0 csin = 0.d0 csinl = 0.d0 cson = 0.d0 csonl = 0.d0 csls = 0.d0 cslo = 0.d0 csom = 0.d0 ampd = 0.d0 avpd = 0.d0 do i = -4,0,1 crdaly(i)=0.d0 end do lqtr = 0 sumcn = 0.d0 sumyn = 0.d0 sumtns = 0.d0 qlno3 = 0.d0 qlnh4 = 0.d0 c*------- calc upper drainage limit (field capacity) water content do i = 1, n udl(i) = thr(i) + (ths(i) - thr(i)) / (1.d0 + 102.d0/at(i)) end do c*------- root length distribution as a function of soil depth c rm = 1.22895d0 c cr = 4.293678d0 c rl = 63.43959d0 c rm = 1.23d0 c cr = 50.0d0 c rl = 56.66d0 c en = dexp(1.d0) do i = 0, n c rdist(i) = rm*(en*z(i)/rl)**cr*dexp(-cr*z(i)/rl) c if (z(i) .le. 35.d0) then c rdist(i) = (-0.0101d0*z(i)+0.4308d0)*z(i)+0.1415d0 c else c rdist(i) = 84.9251d0*dexp(-0.0970144d0*z(i)) c end if rdist(i) = dexp(-0.03d0*z(i)) end do c read in the initial conditions for the water, nh4-n conc., c no3-n conc. and soil water pressure. c if all variables (i.ee. pressure head, nh4-n conc. etc.) are not c available at common soil depths, subroutine idist2 must be used do i=1,7 read(16,*) end do c initial pressure head depths and values call idist2 (p) zrj = p(0) zrp = zrj zr = zrp c*------- set up total head array ht(0) = p(0) do i = 1, n ht(i) = p(i) - z(i) end do ht(np1) = pl - zl p(np1) = pl c initial nh4-n concentration call idist2 (c) c initial no3-n concentration call idist2 (y) c initial soil temperatures call idist2 (t) c*------- set up 15-deg-C lower boundary condition for temperature if (lt .lt. n) then do i = lt, n t(i) = 15.d0 end do end if t(np1) = 15.d0 c*------- calculate nh4 phase exchange coefficient do i = 1, n rkd(i) = rex(i)*rob(i) end do close (15) close (16) print* print* print*,'*********************************************************' print*,'* simulation of nitrogen transformation and transport *' print*,'*********************************************************' print* print* c*------- make initial calls for water properties and write out call wprop do i = 0, np1 thj(i) = th(i) thp(i) = th(i) end do winit = 0.0d0 do i = 1, n winit = winit + th(i) end do winit = winit*dz rday = .true. call output (.true., sofp, soft, wbal) accpln=0.d0 c accevp=0.d0 ar=0.d0 napp=0 napday=.false. klm = 1 maxdwo = iwm allok=.true. c c get 1st day weather, scs curve number, and 1st litter application info c call weather(allok) if (allok) then read(77,*) cn call nextapp(allok) call nextmd endif if (allok) + napday=.not.((kyer(klm).eq.kyr).and.(kdoy(klm).eq.kday)) c********************************************************************** c c ---------------- start the daily simulation loop ----------------- c c********************************************************************** do while ((klm .lt. maxdwo).and.(allok)) do while ((napday).and.(klm.lt.maxdwo)) c* ------ the following code takes the weather data for the day of c* ------ simulation, computes the runoff using the scs curve number c* ----- method, evapotranspiration using hargreves ET. AL, 1980's model iyr = kyer(klm) rain(klm) = crdaly(klm) do na = 1, napp crain(na) = crain(na) + rain(klm) end do c* ---------- compute surface runoff using the scs curve number method c* determine dormant or growing season, set 5 day antecedent rainf. lim if(kdoy(klm) .lt. 90 .or. kdoy(klm) .gt. 305) then rl1 = 1.3d0 rl2 = 2.8d0 else rl1 = 3.6d0 rl2 = 5.3d0 endif c* --------- adjust the curve number for condtions i and iii ---------- if(ar .lt. rl1) then cn1 = 4.2d0*cn/(10.0d0 - 0.058d0*cn) s = 2540.0d0/cn1-25.4d0 elseif(ar .gt. rl2) then cn3 = 23.d0*cn/(10.0d0+0.13d0*cn) s = (2540.0d0/cn3)-25.4d0 else s = (2540.0d0/cn)-25.4d0 endif c* ---------- subtract rainfall greater than 7.62 cm remain = 0.0d0 if(rain(klm) .gt. 7.62d0) then remain = rain(klm) - 7.62d0 rain(klm) = 7.62d0 endif c* -------------------------- compute runoff in cm if(rain(klm) .gt. 0.2d0*s) then runof(klm) = ((rain(klm) - 0.2d0*s)**2)/(rain(klm) + + 0.8d0*s) runof(klm) = runof(klm) + remain else runof(klm) = 0.0d0 endif c* ------- compute rain water and nh4 flux at the soil surface if(rain(klm) .gt. 0.d0) then xinfl(klm) = rain(klm) -runof(klm) +remain else xinfl(klm) = 0.d0 end if cqin = 0.d0 if (xinfl(klm) .gt. 0.d0) then rday = .true. c xtinf = tinf c sflux = xinfl(klm)/xtinf zr = xinfl(klm) zrj = zr zrp = zr p(0) = zr ht(0) = zr c(0) = csin/crdaly(klm) csls = c(0)*xinfl(klm) cslo = c(0)*runof(klm) fcu = csls/xtinf csin = 0.d0 else rday=.false. xinfl(klm) = 0.0d0 c sflux = 0.0d0 c(0) = 0.d0 csls = 0.d0 cslo = 0.d0 fcu = 0.d0 endif c csno3(napp) = 0.0d0 call wprop call etrx c*------- tinfl = time + xtinf c dt = 1.d0 dt = dtt dz = dzz ifc=0 c* ------ if there is no rain, skip this portion and go to statement 26 if (.not. rday) go to 26 c* ---- calculate during rain c if (qu .gt. rx(1)) c + dt = 0.05d0*dz*(ths(1)-th(1))/(qu-rx(1)) c + dt = 0.0005d0*dz*(ths(1)-th(1))/(qu-rx(1)) c + dt = 0.00005d0*dz*(ths(1)-th(1))/(qu-rx(1)) c + dt = 0.000005d0*dz*(ths(1)-th(1))/(qu-rx(1)) c + dt = 0.0000005d0*dz*(ths(1)-th(1))/(qu-rx(1)) c + dt = 0.1d0*dtt c + dt = dtt c if (runof(klm) .gt. 0.d0) dt = 0.05d0*dt c*------- write output files just before rain for redistribution c*------- since last rain c call output (.false., sofp, soft, wbal) do while ((time + 1.1d0*dt) .lt. tinfl) call compute (ifc, noluck, cqin) if (noluck .gt. 0) then write (2, *) 'noluck =', noluck return end if if (dt .gt. dtt) dt = dtt end do if ((tinfl-time) .gt. 0.d0) then dt = tinfl - time call compute (ifc, noluck, cqin) if (noluck .gt. 0) then write (2, *) 'noluck =', noluck return end if time = tinfl end if c*------- write out profile data right after rain c call output (.false., sofp, soft, wbal) c c c ***** cut off infiltration, compute to end of day c do while ((time+0.01d0*dt) .le. timnw) call compute (ifc, noluck, cqin) if (noluck .gt. 0) then write (2, *) 'noluck =', noluck return end if if (dt .gt. dtt) dt = dtt end do if (timnw .gt. (time+1.d-10)) then dt = timnw - time call compute (ifc, noluck, cqin) if (noluck .gt. 0) then write (2, *) 'noluck =', noluck return end if end if go to 40 26 continue c*------- step through no-rain day by dtt do while ((time+0.01d0*dt) .le. timnw) call compute (ifc, noluck, cqin) if (noluck .gt. 0) then write (2, *) 'noluck =', noluck return end if if (dt .gt. dtt) dt = dtt end do if (timnw .gt. (time+1.d-10)) then dt = timnw - time call compute (ifc, noluck, cqin) if (noluck .gt. 0) then write (2, *) 'noluck =', noluck return end if end if 40 continue c*------- write out to files at end of day call output (.true., sofp, soft, wbal) write (*, '(2i4, i6, 2f15.10, f10.5, f13.5, g13.5)') + kyer(klm), kdoy(klm), ifc, dt, res, xinfl(klm), cqin, + wbal/dble(ifc) write (2, '(2i4, i6, 2f15.10, f10.5, f13.5, g13.5)') + kyer(klm), kdoy(klm), ifc, dt, res, xinfl(klm), cqin, + wbal/dble(ifc) write (*,*) dt = dtt time = timnw timnw = timnw + twrite ar=ar+crdaly(klm)-crdaly(klm-5) klm = klm + 1 call weather(allok) if (allok) then napday=.not.((kyer(klm).eq.kyr).and.(kdoy(klm).eq.kday)) else napday=.false. endif end do c*------- end of inner do while ((napday).and.(klm.lt.maxdwo)) if ((klm .ge. maxdwo) .and. (allok)) then c*------- recycle time series arrays c accpln=0.d0 c accevp=0.d0 do ib = -4,0,1 crdaly(ib) = crdaly(klm+(ib-1)) end do kyer(1) = kyer(klm) kdoy(1) = kdoy(klm) tempmx(1) = tempmx(klm) tempmn(1) = tempmn(klm) et(1) = et(klm) pet(1) = pet(klm) runof(1) = runof(klm) xinfl(1) = xinfl(klm) klm = 1 end if if ((allok) .and. .not.(napday)) then c*------- initialize at start of application day napp = napp + 1 c*------- convert from tons/acre litter to ug/cm^2 total N tappn(napp) = percen*aplr*224.265d0 fam(napp) = 0.d0 fav(napp) = 0.d0 crain(napp) = 0.d0 csin = csin + fci*tappn(napp) cson = cson + (1.d0 - fci)*tappn(napp) call nextapp(allok) napday=.true. cav = 0.d0 endif end do c*------- end of outer do while ((klm .lt. maxdwo).and.(allok)) klm = klm - 1 rday = .true. qtr = .true. c myr = kyer(klm) c mdoy = kdoy(klm) call output (.true., sofp, soft, wbal) print *, 'sofp =', sofp, ' soft =', soft c lahey timer routine call timer(iftime) entime = 0.01d0*dble(iftime) c call date_and_time (cp1,cp2,cp3,istr) c entime = 0.001d0*istr(8)+istr(7)+60.d0*(istr(6)+60.d0*(istr(5) c 1 +24.d0*istr(3))) ttime = entime - strt write (*,*)' elapsed execution time= ',ttime/60.,' minutes' print *, char(7) return end c*------- new subroutine amonia c*------- c This function calculates the Arrhenius temperature correction c to a chemical rate involving bacterial action c a = "activation energy" at 25 deg-C (deg-Kelvin) c t = current temperature of the reaction (deg-C) c*------- function arrhen (a, tc) implicit real*8(a-h,o-z) t = tc + 273.15d0 k10 = exp(a*15.d0/(283.15d0*298.15d0)) if (tc .le. 5.d0) then arrhen = 0.d0 else if (tc .lt. 10.d0) then arrhen = k10*(0.2d0*tc - 1.d0) else arrhen = exp(a*(1/t - 1/298.15)) end if return end c*********************************************************************** c this subroutine does nitrogen and water calculations c*********************************************************************** subroutine compute (ifc, noluck, cqin) implicit real*8(a-h,o-z) include 'come4.for' dimension htj(0:izm) c call wprop noluck = 0 ifc = ifc+1 if (ifc .gt. 1.d+5) then noluck = 3 print *, 'ifc exceeds 100,000 in compute' return end if c set up for dt restart do i = 0, n htj(i) = ht(i) end do 10 continue do i = 0, n ht(i) = htj(i) p(i) = ht(i) + z(i) end do call wpset (0,1) cs = 0.5d0*(sc(0)+sc(1)) cb = 0.5d0*(bc(0)+bc(1)) ta = 0.5d0*(at(0)+at(1)) em = 0.5d0*(vm(0)+vm(1)) en = 0.5d0*(vn(0)+vn(1)) cv(i) = dkmh (ck(0), ck(1), p(0), p(1), cs, cb, + ta, em, en, 0.5d0*dz, wv(0)) if (rday .and. (zrj .gt. 0.1d0) .and. (ht(1) .lt. 0.d0)) then alf = 2.d0*cv(0)/dz dtr = dlog((zrj+qe/alf-ht(1))/(qe/alf-ht(1)))/alf if (dt .gt. dtr) dt = dtr end if call pwater (ifc, noluck) if (noluck .gt. 0) return if (noluck .lt. 0) go to 10 call water (ifc, noluck, cqin, dtfl) if (noluck .gt. 0) return if (noluck .lt. 0) go to 10 time = time + dt c call temprt c if (noluck .gt. 0) return c call litter c call amonia c call nitrat c call amonia c*------- reset old water content do i = 1, n thj(i) = th(i) end do c*------- calculate flow to groundwater c gww1 = gww1 + v(n)*dt gww1 = gww1 + (a1*qlp+a2*ql)*dt c*------- adjust dt based on water iterations if (inewt .le. (maxit-2)) then if (p(0) .ge. 0.d0) then dt = 1.062d0*dt else dt = 1.14d0*dt end if if (dt .gt. dtt) dt = dtt else if (inewt .eq. (maxit-1)) then if (p(0) .ge. 0.d0) then dt = 1.001d0*dt else dt = 1.005d0*dt end if if (dt .gt. dtt) dt = dtt else if (inewt .ge. maxit) then dt = 0.619d0*dt end if if (dtfl .gt. 0.d0) then c dt = dtt dt = dmin1(2.d0, dtt) dtfl = 0.d0 end if return end c************************************************************** c* This subroutine calculates root water extraction (cm/h) c* and actual et (cm/h) based on crop stress c************************************************************** subroutine etrx implicit real*8(a-h,o-z) include 'come4.for' c*------- adjust potential et to actual et etr = 0.577d0 evr = 0.0873d0 gk = 12.65d0 c etr = 0.55758d0+4.102d0/gk**1.324d0 c et(klm) = 0.85d0*pet(klm)/24.d0 c etr = 0.686d0 et(klm) = etr*pet(klm)/24.d0 sev = evr*pet(klm)/24.d0 c et(klm) = etr*0.85d0*pet(klm)/24.d0 c sev = etr*0.15d0*pet(klm)/24.d0 if (ck(1) .lt. sev) then qe = ck(1) else qe = sev end if c*------- calculate root extraction, xl = 0.0d0 do i = 1, n rx(i) = rdist(i)*ck(i) xl = xl + rx(i) end do xl = xl*dz c*------- with root stress factor c etk = et(klm) do i = 1, n rst = (th(i) - thr(i)) / (udl(i) - thr(i)) etg = dexp(-dexp(1.932645d0 -gk*rst)) if (p(i) .gt. wilt) then rx(i) = etg*et(klm)*rx(i)/xl else rx(i) = 0.0d0 end if end do return end c********************************************************************** c* this subroutine calculates both actual and potential c* evapotranspiration c********************************************************************** subroutine evapo implicit real*8(a-h,o-z) include 'come4.for' c* ---------- convert temperatures to degrees f ----------------- tavef = 0.5d0*(tmaxf + tminf) c* ------------ compute potential evapotranspiration ------------- eo1 = 0.0075 * rad*tavef*xkt*((tmaxf-tminf)**(.5)) c* ----------- convert langleys/day to mm/day ------------------- c* ------------ and adjust for temperature ---------------------- eeq = eo1 / 58.5d0 eeq = eeq/10.d0 pet(klm) = eeq*1.1d0 if (tmax .gt. 35.d0) + pet(klm) = eeq*((tmax-35.d0)*0.05d0+1.1d0) if (tmax .lt. 5.d0) + pet(klm) = eeq*0.01d0*dexp(.18d0*(tmax+20.d0)) c* --------- compute actual evapotranspiration ----------------- sev = 0.15d0*pet(klm)/24.d0 return end c*********************************************************************** c* this subroutine sets up the initial soil water pressure, soil water * c* content, nh4-n content, and no3-n content. linear interpolation was * c* used to predict the value of these variables at any soil depth from * c* input values at specified soil depths. * c*********************************************************************** subroutine idist1 (sval) implicit real*8(a-h,o-z) include 'come4.for' real*8 sval(0:izm), xxx(0:izm), c1(0:izm) c read number of data points read(15, '(i3)') nin c read soil depth locations do i = 1, np1 xxx(i) = 0.d0 end do read(15,*) (xxx(i),i=1,nin) c read data values for corresponding depths read(15,*) (c1(i),i=1,nin) i = 1 do k = 0, np1 a = z(k) 5 if(a .le. xxx(i+1)) then sval(k)=c1(i)+(a-xxx(i))*((c1(i+1)-c1(i))/(xxx(i+1)-xxx(i))) goto 10 else if (xxx(i) .gt. xxx(i+1)) then sval(k) = c1(i)+(a-xxx(i))*((c1(i)-c1(i-1)) + /(xxx(i)-xxx(i-1))) goto 10 end if i = i + 1 goto 5 10 continue end do return end c*********************************************************************** c* this subroutine sets up the initial soil water pressure, soil water * c* content, nh4-n content, and no3-n content. linear interpolation was * c* used to predict the value of these variables at any soil depth from * c* input values at specified soil depths. * c*********************************************************************** subroutine idist2 (sval) implicit real*8(a-h,o-z) include 'come4.for' real*8 sval(0:izm), xxx(0:izm), c1(0:izm) c read number of data points read(16,*) nin c read soil depth locations do i = 1, np1 xxx(i) = 0.d0 end do read(16,*) (xxx(i),i=1,nin) c read data values for corresponding depths read(16,*) (c1(i),i=1,nin) i = 1 do k = 0, np1 a = z(k) 5 if(a .le. xxx(i+1)) then sval(k)=c1(i)+(a-xxx(i))*((c1(i+1)-c1(i))/(xxx(i+1)-xxx(i))) goto 10 else if (xxx(i) .gt. xxx(i+1)) then sval(k) = c1(i)+(a-xxx(i))*((c1(i)-c1(i-1)) + /(xxx(i)-xxx(i-1))) goto 10 end if i = i + 1 goto 5 10 continue end do return end c********************************************************************** c* c* Given z(i) and variable v(i), interpolate the value, val, at zv c* c********************************************************************** subroutine interp (z, v, zv, val) implicit real*8(a-h,o-z) parameter (izm=107) double precision z(0:izm), v(0:izm) common /pr4/ lt, nm1, n, np1, inewt, maxit do i = 0, n if (zv .le. z(i+1)) then val = v(i) +(zv-z(i))*(v(i+1)-v(i))/(z(i+1)-z(i)) go to 10 end if end do val = v(n) + (zv-z(n))*(v(n)-v(nm1))/(z(n)-z(nm1)) 10 continue return end c********************************************************************** c* this subroutine keeps track of the day of the year, thus if the c* dataset includes the weather data and the corresponding day of the c* year this subroutine is skipped. c********************************************************************** subroutine jconvt implicit real*8(a-h,o-z) include 'come4.for' dimension ic(12) data ic/31,28,31,30,31,30,31,31,30,31,30,31/ quo = 0.25d0*kyer(klm) ipt = idint(quo) rplt = float(ipt) rem = quo - rplt if(rem .ne. 0.0) then ic(2)=28 else ic(2)=29 endif jul = kdoy(klm) i = 1 10 if(jul - ic(i) .le. 0) goto 20 jul = jul - ic(i) i = i + 1 goto 10 20 imm = i idd = jul return end c********************************************************************** c* this subroutine calculates the mineralization and volatilization of c* the litter nitrogen, and its conversion to soil OM (organic matter) c* c********************************************************************** subroutine litter implicit real*8(a-h,o-z) include 'come4.for' real*8 rk(3) c +, ak data rk /0.049d0,0.019d0,0.0060d0/ c +, ak /-6750.d0/ tc = blt do i = 1, napp if (fam(i) .lt. 0.85d0) then c*------- the mineralization process rate constant fonr = 1.d0 - fam(i) if ((1.d0-fam(i)) .gt. 0.8d0) then akm = rk(1) else if ((1.d0-fam(i)) .gt. 0.65d0) then akm = rk(2) else akm = rk(3) end if akm = akm/24.d0 c*------- water content correction if (th(0) .lt. 0.1d0) then ae=9.d0 af=0.d0 elseif (th(0) .lt. 0.2d0) then ae=1.d0 af=0.8d0 else if (th(0) .lt. 0.668d0) then ae=-2.13675d0 af=1.42735d0 else ae = 0.d0 af = 0.d0 endif akm=akm*(ae*th(0)+af) c*------- air temp adjustment to bacterial action akm = tc*akm c*------- mineralization dfam = fonr*(1.d0 - dexp(-akm*dt)) if ((fam(i)+dfam) .ge. 0.85d0) dfam = 0.85d0 - fam(i) fam(i) = fam(i) + dfam d3 = dfam*tappn(i) csin = csin + d3 ampd = ampd + d3 cson = cson - d3 c*------- calculate volatilization if any inorganic-N left dfav = 0.d0 if (fav(i) .lt. fam(i)) then fmaxv = 0.01d0*( 5.69783d0 +0.229995d0*tave + +0.000459594d0*tappn(i) -0.353627d0*crain(i) ) akv = -0.0170028d0 +9.29487d-4*tave -2.41901d-7*tappn(i) + +0.00108953d0*crain(i) if ((akv .gt. 0.d0) .and. (fmaxv .gt. 0.d0)) then akv = tc*akv fmaxv = tc*fmaxv dfav = (fmaxv - fav(i))*(1.d0 - dexp(-akv*dt)) if (dfav .lt. 0.d0) dfav = 0.d0 if ((fav(i)+dfav) .ge. fam(i)) dfav = fam(i) -fav(i) fav(i) = fav(i) + dfav end if end if d3 = dfav*tappn(i) c if (d3 .gt. csin) d3 = csin avpd = avpd + d3 cav = cav + d3 csin = csin - d3 else go to 10 end if c*------- if only 15% of total applied N left, dump remainder to soil c* organic matter, and reset litter to zero if (fam(i) .ge. 0.85d0) then don = (1.d0 - fam(i))*tappn(i) cson = cson - don csom = csom + don end if end do 10 continue cslr = csin + cson return end c*********************************************************************** c c This subroutine reads information regarding the next poultry litter c application and verifies compatibility with the current weather data c c*********************************************************************** subroutine nextapp (allok) implicit real*8(a-h,o-z) include 'come4.for' logical allok read(77,*,end=89) kyr, kday, aplr, percen, fci if (kyr.lt.kyer(klm)) then if (.not.((kyr.eq.0).and.(kyer(klm).eq.99))) then write (*,*)' Litter Application Date preceeds Weather Date' write (*,*)' Simulation Ended' allok=.false. endif elseif (kyr.eq.kyer(klm)) then if (kday.lt.kdoy(klm)) then write (*,*)' Litter Application Date preceeds Weather Date' write (*,*)' Simulation Ended' allok=.false. endif endif return 89 kday=375 end c*********************************************************************** c c This subroutine reads information regarding the next c measurement day for which model output will be generated c c*********************************************************************** subroutine nextmd implicit real*8(a-h,o-z) include 'come4.for' read(1,*,end=89) myr, mdoy, (pm(i), i=1,5), tm(1), tm(2), wm(1), + wm(2) 89 continue return end c*------- new subroutine nitrat c**************************************************************************** c* this subroutine provides the solution to the nitrate transport and * c* transformation equation under transient-flow conditions. it also * c* calculates the nitrate uptake by plant roots. the method of solution * c* is the finite difference approximation method. the rate of nitrification * c* denitrification and the distribution of ammonium exchange are obtained * c* from functions, zz1, zz2, and zz3, respectively. * c**************************************************************************** c subroutine nitrat c**************************************************************************** c* this subroutine prints the results from the model predictions for soil * c* water, nh4-n, and no3-n at specified times. it also calculates the total * c* amounts of nitrogen in the soil system and that taken up by plants. * c**************************************************************************** subroutine output (eod, sofp, soft, wbal) implicit real*8(a-h,o-z) include 'come4.for' logical eod c* ------ write out the predicted soil water pressure, soil water c* ------ content, flow rate, nh4-n conc. and no3-n conc. at each c* ----------------------- soil depth c*------- 12 --> *1.out c if (rday) then c if ((kyer(klm) .eq. myr) .and. (kdoy(klm) .eq. mdoy) .and. c + eod) then if ((time .gt. 0.d0) .and. eod .and. qtr) then timm = time/24.0d0 c write (12, *) char(12) c write (12,200) timm c write (12,199) kyer(klm),imm,idd c write (12,299) c write (12,301) 199 format('1',30x,'yy-mm-dd =', i3,i3,i4/) 200 format('1',30x,'time, days =',f12.4/) c 299 format( c +' depth pre. head soil water flow vel nh4 conc no3 co c +nc temp cond') c 301 format( c +' cm cm cm3/cm3 cm/day ug/cm^3 ug/cm c +^3 c cm/day'/) 299 format( +' depth pre. head soil water flow vel temp cond') 301 format( +' cm cm cm3/cm3 cm/day c cm/day'/) wflx = 0.0d0 zp = 0.d0 do 20 i = 0 ,n if (z(i) .ge. zp) then c*------- note: qu and wflx are not fully determined; need to be fixed if (i .eq. 0) then c wflx = qu wflx = qin else wflx = -cv(i-1)*(ht(i)-ht(i-1))/dz end if call interp (z, p, zp, pp) call interp (z, th, zp, thp) call interp (z, c, zp, cp) call interp (z, y, zp, yp) call interp (z, t, zp, tp) call interp (z, ck, zp, ckp) c*------- 12 --> *1.out c write(12,499) zp, pp, thp, wflx*24, cp, yp, tp, ckp*24 c write (12,499) zp, pp, thp, wflx*24, tp zp = zp + dzp end if 499 format(f6.1, 2x, f10.2, 5x, f8.4, 5x, f8.4, 4x, f8.3, 2x, f8.3, + 2x, f6.2, 3f10.5) 20 continue wflx = -ck(n)*(ht(np1)-ht(n))/dzl c write(12,499) z(np1), p(np1), th(np1), wflx*24, c(np1), c + y(np1), t(np1), ck(np1)*24 c write (12,499) z(np1), p(np1), th(np1), wflx*24, c + t(np1) c WRITE (12, *) ' RES =', RES, ' DT =', DT, ' V(N) =', V(N), c + ' PLNW =', PLNW end if c* ------- calculate by integrating over the whole soil profile, the c* ------- nh4-n conc. in solution, in exchangeable phase and total n. c* ------------- also, the total no3-n in the profile is given. c do i = 1, n c aa(i) = th(i)*c(i) c bb(i) = th(i)*y(i) c cc(i) = rob(i)*rex(i) c end do c call qsf(dz, aa, xnh4, n) c call qsf(dz, bb, xno3, n) c call qsf(dz, cc, enh4, n) 30 continue c* -------- this logic computes the mass balance for soil water wfin = 0.0d0 do i = 1, n wfin = wfin + th(i) end do wfin = wfin*dz plnup = plnw-plnw1 c wevap = pet(klm) c accevp = accevp + wevap dev = cev - cev1 gwater = gww1 - gww2 wbal = winit -wfin +crdaly(klm) -runof(klm) -dev -gwater -plnup c* ------- this logic computes the mass balance for nitrogen c xtnh4 = xnh4 + enh4 c tnfin = xno3 + xtnh4 c gwn1 = gwno3 + gwnh4 c gwn = gwn1 - gwn2 c pln = plno3 + plnh4 - pln1 c dnit = dnitrf - dnitr1 c tnbal = tninit - tnfin + csls - pln - dnit - gwn c*------- end of day additions to plant N and total N in solution sums c if (eod) then c accpln = accpln + pln c sumtns = sumtns + xno3 + xnh4 c end if c*------- 12 --> *1.out c* - write the mass balance for soil water and nitrogen to output files c*------- 14 --> *2.out soil profile nitrogen data c*------- 25 --> *3.out soil profile water data if ((time .gt. 0.0d0) .and. eod) then write (25,23) iyr,imm,idd,nday,winit,wfin,crdaly(klm), + xinfl(klm),dev,runof(klm),gwater,plnup, pet(klm), wbal end if 23 format(1x, i2, 2i3, i4, 9f9.3, g13.3) c*------- 18 --> *4.out litter nitrogen changes if ((time .gt. 0.d0) .and. eod) then c write (18, 24) iyr, nday, napp, csonl/10.d0, cson/10.d0, c + ampd/10.d0, csinl/10.d0, csin/10.d0, csls/10.d0, cslo/10.d0, c + c(0), avpd/10.d0, cav/10.d0, csom/10.d0, tave end if 24 format (i2, 2i4, 12f9.3) c*------- if it's time, calculate quarterly data output c*------- 7 --> *5.out M. Cochran's data C*------- write out modeled sensor readings on meas. day 2 --> *6.out if ((kyer(klm) .eq. myr) .and. (kdoy(klm) .eq. mdoy) .and. + (time .ne. 0.d0) .and. eod) then call interp (z, p, 30.d0, p1) call interp (z, p, 60.d0, p2) call interp (z, p, 90.d0, p3) call interp (z, p, 120.d0, p4) call interp (z, p, 200.d0, p5) call interp (z, t, 10.16d0, t1) call interp (z, t, 30.48d0, t2) call interp (z, th, 10.16d0, th1) call interp (z, th, 30.48d0, th2) c write (2, 999) myr, mdoy, p1, p2, p3, p4, p5, t1, t2, th1, th2 999 format(1x, i2, i4, 5f9.2, 2f8.2, 2f8.5) c*------- write out objective functions for measurement day, c* 3 --> *7.out ofp(1) = p1-pm(1) ofp(2) = p2-pm(2) ofp(3) = p3-pm(3) ofp(4) = p4-pm(4) ofp(5) = p5-pm(5) ofp(6) = p1-p2-(pm(1)-pm(2)) ofp(7) = p2-p3-(pm(2)-pm(3)) ofp(8) = p3-p4-(pm(3)-pm(4)) ofp(9) = p4-p5-(pm(4)-pm(5)) ofp(10) = 1000.d0*(th1 - wm(1)) ofp(11) = 1000.d0*(th2 - wm(2)) c sofp = 0.d0 do i = 1, 11 sofp = sofp + ofp(i)*ofp(i) end do c write (3, '(i3, i4, 12e18.7)') myr, mdoy, (ofp(i),i=1,11), sofp print *, myr, mdoy, (ofp(i), i=1,11), sofp oft(1) = t1-tm(1) oft(2) = t2-tm(2) oft(3) = t1-t2-(tm(1)-tm(2)) soft = soft + oft(1)*oft(1) + oft(2)*oft(2) + oft(3)*oft(3) soft = soft + 0.25d0*(t(n)-15.d0)**2 c write (4, '(i3, i4, 4e18.7)') myr, mdoy, oft(1), oft(2), c + oft(3), soft print *, myr, mdoy, oft(1), oft(2), oft(3), soft c print *, call nextmd end if c*------- if not end of day, then don't reset variables if (.not. eod) return c* --------- reinitialize the variables for the new day simulation pln1 = plno3 + plnh4 plno31 = plno3 plnh41 = plnh4 plnw1 = plnw cev1 = cev dnitr1 = dnitrf tninit = tnfin winit = wfin gwn2 = gwn1 gww2 = gww1 cavl = cav csinl = csin csonl = cson ampd = 0.d0 avpd = 0.d0 c if (time .gt. 0.d0) write (*,*) kyer(klm), kdoy(klm), ifc, dt c write (*,*) return end c********************************************************************** c* This integration subroutine uses simple summation*dz c* c********************************************************************** c subroutine qsf (dz,y,zint,n) c*********************************************************************** c* * c* this routine reads in the weather data from the dataset * C* * c*********************************************************************** subroutine readit (allok) implicit real*8(a-h,o-z) include 'come4.for' logical allok c* -------- read in weather data which should be in the formatt ------ c* -------- described in the supporting document. --------------------- c crdaly(klm)= 0.0d0 read(20,1052,end=121,err=121) kyer(klm),kdoy(klm),tmaxf, + tminf, crdaly(klm) 1052 format(i2,i3,f3.0,f3.0,f5.0) c*------ kill program if extra line past weather data read, or yr 2000 c*------ needs to be fixed in four years c if (kyer(klm) .lt. kyer(klm-1)) go to 121 c* --------- change temperatures read from file from degrees f to ---- c* --------- degrees c tmax = (tmaxf - 32.0d0) * 5.0d0/9.0d0 tmin = (tminf - 32.0d0) * 5.0d0/9.0d0 tempmx(klm) = tmax tempmn(klm) = tmin c* ------------ conversion of daily precip. to cm crdaly(klm) = crdaly(klm)*0.00254d0 if(crdaly(klm) .le. 0.254d0) crdaly(klm) = 0.0d0 return 121 continue write (*,*) ' End of Weather Data - Simulation Ended' allok=.false. return end c*------- new subroutine temprt C********************************************************************** c* this subroutine provides the solution to the transient heat flow c* equation. the thermal diffusivity was computed as function of soil c* water content. c********************************************************************** subroutine temprt implicit real*8(a-h,o-z) include 'come4.for' c rb = dt/dz c ra = rb/dz c rd = 2.d0*dt/(dz*dz*3.d0) c*------- set upper boundary condition t(0) = tave dzm = 0.5d0*(dzg + 1.5d0*dz) rd = dt/dzm dm = rd/(dzm/tdg + dz/tdif(1)) aa(1) = 0.d0 cc(1) = -0.5d0*rd*(tdif(1)+tdif(2))/dz bb(1) = 1.d0 -cc(1) +dm r(1) = t(1) + dm*tave c*------- calculate difference matrix, Euler implicit time step rd = 0.5d0*dt/(dz*dz) tdi = tdif(1) tdp = tdif(2) c do i = 2, nm1 do i = 2, lt-1 tdm = tdi tdi = tdp tdp = tdif(i+1) aa(i) = -rd*(tdm+tdi) cc(i) = -rd*(tdi+tdp) bb(i) = 1.d0 -aa(i) -cc(i) r(i) = t(i) end do c*------- set lower boundary condition c t(np1) = 15.d0 if (lt .le. n) then rd = 0.5d0*dt/dz dm = rd*(tdif(lt-1) + tdif(lt))/dz aa(lt-1) = -rd*(tdif(lt-2) + tdif(lt-1))/dz cc(lt-1) = 0.d0 bb(lt-1) = 1.d0 -aa(lt-1) +dm r(lt-1) = t(lt-1) +tl*dm else rd = dt/dz dm = rd*tdif(n)/(ztl-cl) aa(n) = -0.5d0*rd*(tdif(n-1) + tdif(n))/dz cc(n) = 0.d0 bb(n) = 1.d0 -aa(n) +dm r(n) = t(n) +tl*dm end if c*------- solve difference method call dgtsl (lt-1, aa, bb, cc, r, info) if (info .gt. 0) then noluck = 2 return end if do i = 1, lt-1 t(i) = r(i) end do c*------- modifiy Imax in Michaelis relation by temperature do i = 1, n if (t(i) .le. 10.d0) then qme(i) = 0.d0 else if ((t(i) .gt. 10.d0) .and. (t(i) .lt. 15.d0)) then qme(i) = qm*(0.2d0*(t(i) - 10.d0)) else qme(i) = qm end if c qme(i) = 0.d0 end do return end subroutine dgtsl(n,c,d,e,b,info) implicit double precision (a-h,o-z) parameter (izm=107) integer n,info double precision c(0:izm), d(0:izm), e(0:izm), b(0:izm) c c dgtsl given a general tridiagonal matrix and a right hand c side will find the solution. c c on entry c c n integer c is the order of the tridiagonal matrix. c c c double precision(n) c is the subdiagonal of the tridiagonal matrix. c c(2) through c(n) should contain the subdiagonal. c on output c is destroyed. c c d double precision(n) c is the diagonal of the tridiagonal matrix. c on output d is destroyed. c c e double precision(n) c is the superdiagonal of the tridiagonal matrix. c e(1) through e(n-1) should contain the superdiagonal. c on output e is destroyed. c c b double precision(n) c is the right hand side vector. c c on return c c b is the solution vector. c c info integer c = 0 normal value. c = k if the k-th element of the diagonal becomes c exactly zero. the subroutine returns when c this is detected. c c linpack. this version dated 08/14/78 . c jack dongarra, argonne national laboratory. c c no externals c fortran dabs c c internal variables c integer k,kb,kp1,nm1,nm2 double precision t c begin block permitting ...exits to 100 c info = 0 c(1) = d(1) nm1 = n - 1 if (nm1 .lt. 1) go to 40 d(1) = e(1) e(1) = 0.0d0 e(n) = 0.0d0 c do 30 k = 1, nm1 kp1 = k + 1 c c find the largest of the two rows c if (dabs(c(kp1)) .lt. dabs(c(k))) go to 10 c c interchange row c t = c(kp1) c(kp1) = c(k) c(k) = t t = d(kp1) d(kp1) = d(k) d(k) = t t = e(kp1) e(kp1) = e(k) e(k) = t t = b(kp1) b(kp1) = b(k) b(k) = t 10 continue c c zero elements c if (c(k) .ne. 0.0d0) go to 20 info = k c ............exit go to 100 20 continue t = -c(kp1)/c(k) c(kp1) = d(kp1) + t*d(k) d(kp1) = e(kp1) + t*e(k) e(kp1) = 0.0d0 b(kp1) = b(kp1) + t*b(k) 30 continue 40 continue if (c(n) .ne. 0.0d0) go to 50 info = n go to 90 50 continue c c back solve c nm2 = n - 2 b(n) = b(n)/c(n) if (n .eq. 1) go to 80 b(nm1) = (b(nm1) - d(nm1)*b(n))/c(nm1) if (nm2 .lt. 1) go to 70 do 60 kb = 1, nm2 k = nm2 - kb + 1 b(k) = (b(k) - d(k)*b(k+1) - e(k)*b(k+2))/c(k) 60 continue 70 continue 80 continue 90 continue 100 continue c return end c********************************************************************** c* this subroutine provides the solution to the tridiagonal matrix- c* vector equations for subroutines water, amonia and nitrat using c* thomas algorithm. c********************************************************************** subroutine tridia (nn, a, b, c, d) implicit double precision (a-h,o-z) parameter (izm=107) dimension a(0:izm), b(0:izm), c(0:izm), d(0:izm) do n = 2 , nn r = a(n) / b(n - 1) b(n) = b(n) - r * c(n - 1) d(n) = d(n) - r * d(n - 1) end do d(nn) = d(nn) / b(nn) do n = nn - 1 , 1 , -1 d(n) = (d(n) - c(n) * d(n + 1)) / b(n) end do return end c******************************************************************************* c* this subroutine provides the numerical solution to the water flow equation. * c* the method of solution is based on a finite difference approximation c******************************************************************************* subroutine pwater (ifc, noluck) implicit real*8(a-h,o-z) include 'come4.for' dimension hti(0:izm) c* --------- setup the tridiagonal matrix of the system of linear c* --------- equations resulting from finite difference approximation c* --------- to the water flow equation using an Euler time step c*------- calculate et, root extraction and surface evaporation noluck = 0 pdt = prk*dt if (rday .and. (zrj .gt. 0.d0)) then ur = 0.5d0 urmax = 1.5d0 else ur = 0.8d0 urmax = 3.d0 end if url = ur call wprop call etrx do i = 0, np1 thp(i) = th(i) hti(i) = ht(i) end do c*------- iterate by Newton's method to solve Richards' equation do 60 inewt = 1, maxit c*------- set up tridiagonal solution to Richards' equation pdt = prk*dt c pdt = dt rb = pdt/dz ra = rb/dz c*------- restart point for iteration 10 continue c*------- upper boundary condition aa(1) = 0.0d0 cc(1) = -ra*cv(1) bb(1) = ch(1) - cc(1) if (zrj .le. 0.d0) then qinp = -qe zrp = ht(1) +0.5d0*qe*dz/cv(0) ri = thp(1)-thj(1) -ra*(cv(1)*(ht(2)-ht(1))) + +rb*qe+rx(1)*pdt else zrp = (zrj -qe*pdt +2.d0*rb*cv(0)*ht(1)) + /(1.d0+2.d0*rb*cv(0)) zr = ((zrj-qe*dt)*dz +2.d0*dt*cv(0)*ht(1)) + / (dz +2.d0*dt*cv(0)) if ((zrp .lt. 0.d0) .or. (zr .lt. 0.d0)) then c if (zrp .lt. 0.d0) then c zrp = 0.d0 qinp = zrj/dt -qe zrp = ht(1) +0.5d0*dz*qinp/cv(0) else cc(1) = -4.d0*ra*cv(1)/3.d0 bb(1) = ch(1) - cc(1) +8.d0*ra*cv(0)/3.d0 c qinp = 2.d0*cv(0)*(zrj-qe*pdt-ht(1))/(dz+2.d0*pdt*cv(0)) qinp = -2.d0*cv(0)*(ht(1)-zrp)/dz end if ri = thp(1)-thj(1) -ra*(cv(1)*(ht(2)-ht(1))) + +rb*(qe-qinp)+rx(1)*pdt end if ht(0) = zrp p(0) = zrp r(1) = -ri call wpset (0, 0) c*------- Jacobian enhancements c cc(1) = cc(1) +0.25d0*rb*(dk(1)+dk(2)) c bb(1) = bb(1) +0.25d0*rb*(dk(2)) cki = ck(1) ckp = ck(2) wv2 = wv(1) dh2 = ht(2) - ht(1) dp2 = p(2) - p(1) cv2 = cv(1) dki = dk(0) dkp = dk(1) do 20 i = 2 ,nm1 dkm = dki dki = dkp dkp = dk(i+1) cv1 = cv2 cv2 = cv(i) dp1 = dp2 dp2 = p(i+1) - p(i) dh1 = dh2 dh2 = ht(i+1) - ht(i) wv1 = wv2 wv2 = wv(i) ckm = cki cki = ckp ckp = ck(i+1) ri = thp(i)-thj(i) -ra*(cv2*dh2 - cv1*dh1) +rx(i)*pdt if (dabs(ri) .lt. 1.d-200) ri = dsign(1.d-200,ri) r(i) = -ri aa(i) = -ra*cv1 bb(i) = ch(i) + ra*(cv1+cv2) cc(i) = -ra*cv2 c*------- Jacobian enhancements c aa(i) = aa(i) -rb*0.25d0*(dkm+dki) c bb(i) = bb(i) +rb*0.25d0*(dkp-dkm) c cc(i) = cc(i) +rb*0.25d0*(dki+dkp) 20 continue c*------- lower boundary condition ri = thp(n)-thj(n) -rb*ck(n)*(ht(np1)-ht(n))/dzl + +ra*cv(nm1)*(ht(n)-ht(nm1)) +rx(n)*pdt r(n) = -ri aa(n) = -ra*cv(nm1) bb(n) = ch(n) -aa(n) +rb*ck(n)/dzl cc(n) = 0.0d0 c*------- Jacobian enhancements c aa(n) = aa(n) -rb*0.25d0*(dk(nm1)+dk(n)) c bb(n) = bb(n) +rb*(dk(n) -0.25d0*(dk(nm1)+dk(n))) c*------- test well-posedness of tridiagonal solution c do i = 1, n c if (dabs(bb(i)) .eq. 0.d0) noluck = 2 c if (dabs(bb(i)) .lt. (dabs(aa(i)) + dabs(cc(i))) ) c + noluck = 2 c if (noluck .eq. 2) then c print *, 'poorly posed tridiag in water at i =', i c return c end if c end do c*------- ensure diagonal dominance of tridiagonal solution c do i = 1, n c dj = dabs(aa(i)) + dabs(cc(i)) - dabs(bb(i)) c if (dj .gt. 0.d0) then c bb(i) = bb(i) + dsign(dj, bb(i)) c end if c end do c*-------- call thomas algorithm to solve for soil water pressure call dgtsl (n, aa, bb, cc, r, info) if (info .gt. 0) then noluck = 2 return end if c*------- advance p(i) by Newton correction do i = 1, n c p(i) = p(i) + r(i) p(i) = p(i) + ur*r(i) c if(p(i) .gt. pmax) p(i) = pmax ht(i) = p(i) - z(i) end do if (zrj .le. 0.d0) then zrp = ht(1) +0.5d0*qe*dz/cv(0) else zrp = (zrj -qe*pdt +2.d0*rb*cv(0)*ht(1)) + /(1.d0+2.d0*rb*cv(0)) zr = ((zrj-qe*dt)*dz +2.d0*dt*cv(0)*ht(1)) + / (dz +2.d0*dt*cv(0)) if ((zrp .lt. 0.d0) .or. (zr .lt. 0.d0)) then qinp = zrj/dt -qe zrp = ht(1) +0.5d0*qinp*dz/cv(0) end if end if ht(0) = zrp p(0) = zrp c* --------- update soil water properties call wprop c*------- calculate et, root extraction and surface evaporation call etrx do i = 0, n thp(i) = th(i) end do c*------- calculate residual mass balance (cm^3/cm^2) if (zrj .le. 0.d0) then qinp = -qe zrp = ht(1) - 0.5d0*qe*dz/cv(0) else zrp = (zrj -qe*pdt +2.d0*rb*cv(0)*ht(1)) + /(1.d0+2.d0*rb*cv(0)) zr = ((zrj-qe*dt)*dz +2.d0*dt*cv(0)*ht(1)) + / (dz +2.d0*dt*cv(0)) if ((zrp .lt. 0.d0) .or. (zr .lt. 0.d0)) then c if (zrp .lt. 0.d0) then c zrp = 0.d0 qinp = zrj/dt -qe zrp = ht(1) +0.5d0*dz*qinp/cv(0) else c qinp = 2.d0*cv(0)*(zrj-qe*pdt-ht(1))/(dz+2.d0*pdt*cv(0)) qinp = -2.d0*cv(0)*(ht(1)-zrp)/dz end if end if ht(0) = zrp p(0) = zrp ql = -ck(n)*(ht(np1)-ht(n))/dzl qlp = ql qnp = qinp - qlp sdth = 0.d0 sresp = 0.d0 do i = 1, n sdth = sdth + thp(i)-thj(i) qnp = qnp - rx(i)*dz sresp = sresp + dabs(r(i)*1.d-2) end do res1 = sdth*dz - qnp*pdt res2 = sresp resp = dsqrt(dabs(res1)*res2) c resp = res1 if (inewt .eq. 1) then respl = resp url = ur else if (dabs(resp) .le. dabs(respl)) then ur = ur + 0.1d0*(urmax-ur) do i = 0, n hti(i) = ht(i) end do else if (ur .gt. 0.1d0) then ur = dmax1(0.1d0, 0.5d0*ur) do i = 0, n ht(i) = hti(i) p(i) = ht(i) + z(i) end do call wprop do i = 0, n thp(i) = th(i) end do call etrx go to 10 end if end if write (*, '(1h+, i7, a, i7, 2f15.10, f6.3)') + inewt, ' NEWT IFC', ifc, resp, dt, ur c if (dt .ge. 0.1d0) then if (dabs(resp) .lt. 1.d-6) go to 70 c else c if (dabs(resp) .lt. 1.d-5) go to 70 c end if respl = resp 60 continue 70 continue if (dt .lt. 1.d-100) then noluck = 4 print *, 'dt goes to zero in subroutine pwater' return end if if (inewt .gt. maxit) then dt = 0.619*dt noluck = -1 if (dt .lt. 1.d-100) then print *, 'endless loop in subroutine pwater' noluck = 1 return end if return end if c* --------- calculate water extracted by roots from the root zone srxp = 0.d0 do i = 1, n srxp = srxp + rx(i) end do return end c******************************************************************************* c* this subroutine provides the numerical solution to the water flow equation. * c* the method of solution is based on a finite difference approximation c******************************************************************************* subroutine water (ifc, noluck, cqin, dtfl) implicit real*8(a-h,o-z) include 'come4.for' dimension hti(0:izm) c* --------- setup the tridiagonal matrix of the system of linear c* --------- equations resulting from finite difference approximation c* --------- to the water flow equation using an Euler time step c*------- calculate et, root extraction and surface evaporation noluck = 0 dtfl = 0 if (rday .and. (zrp .gt. 0.d0)) then ur = 0.5d0 urmax = 1.5d0 else ur = 0.8d0 urmax = 3.d0 end if url = ur call wprop call etrx do i = 0, np1 hti(i) = ht(i) end do c*------- iterate by Newton's method to solve Richards' equation do 60 inewt = 1, maxit c*------- set up tridiagonal solution to Richards' equation rb = a2*dt/dz ra = rb/dz ap = a1/prk c ap = 0.d0 c*------- restart point for iteration 10 continue c*------- upper boundary condition aa(1) = 0.0d0 cc(1) = -ra*cv(1) bb(1) = ch(1) - cc(1) c if (zrp .le. 0.d0) then if (zrj .le. 0.d0) then qin = -qe zr = ht(1) - 0.5d0*qe*dz/cv(0) ri = th(1)-thj(1) -ap*(thp(1)-thj(1)) + -ra*cv(1)*(ht(2)-ht(1)) +rb*qe + +a2*rx(1)*dt else c zr = ( zrj +ap*(zrp-zrj) -a2*qe*dt +2.d0*rb*cv(0)*ht(1) ) c + /(1.d0+2.d0*rb*cv(0)) c if (zr .lt. 0.d0) then c zr = 0.d0 c qin = zrp/dt -qe c else c qin = -2.d0*cv(0)*(ht(1)-zr)/dz c end if if (qinp .eq. zrj/dt-qe) then zr = ht(1) +0.5d0*dz*qinp/cv(0) qin = qinp else cc(1) = -4.d0*ra*cv(1)/3.d0 bb(1) = ch(1) - cc(1) +8.d0*ra*cv(0)/3.d0 zr = ( zrj +ap*(zrp-zrj) -a2*qe*dt +2.d0*rb*cv(0)*ht(1) ) + /(1.d0+2.d0*rb*cv(0)) qin = -2.d0*cv(0)*(ht(1)-zr)/dz end if ri = th(1)-thj(1) -ap*(thp(1)-thj(1)) + -ra*cv(1)*(ht(2)-ht(1)) +rb*(qe-qin) + +a2*rx(1)*dt end if ht(0) = zr p(0) = zr r(1) = -ri call wpset (0, 0) c*------- Jacobian enhancements c cc(1) = cc(1) +0.25d0*rb*(dk(1)+dk(2)) c bb(1) = bb(1) +0.25d0*rb*dk(2) cki = ck(1) ckp = ck(2) wv2 = wv(1) dh2 = ht(2) - ht(1) dp2 = p(2) - p(1) dki = dk(0) dkp = dk(1) cv2 = cv(1) do 20 i = 2 ,nm1 dkm = dki dki = dkp dkp = dk(i+1) dp1 = dp2 dp2 = p(i+1) - p(i) dh1 = dh2 dh2 = ht(i+1) - ht(i) wv1 = wv2 wv2 = wv(i) ckm = cki cki = ckp ckp = ck(i+1) cv1 = cv2 cv2 = cv(i) ri = th(i)-thj(i) -ap*(thp(i)-thj(i)) + -ra*( cv2*dh2 - cv1*dh1 ) +a2*rx(i)*dt if (dabs(ri) .lt. 1.d-200) ri = dsign(1.d-200,ri) r(i) = -ri aa(i) = -ra*cv1 cc(i) = -ra*cv2 bb(i) = ch(i) - aa(i) -cc(i) c*------- Jacobian enhancements c aa(i) = aa(i) -rb*0.25d0*(dkm+dki) c bb(i) = bb(i) +rb*0.25d0*(dkp-dkm) c cc(i) = cc(i) +rb*0.25d0*(dki+dkp) 20 continue c*------- lower boundary condition ri = th(n)-thj(n) -ap*(thp(n)-thj(n)) + -rb*ck(n)*(ht(np1)-ht(n))/dzl + +ra*cv(nm1)*(ht(n)-ht(nm1)) +a2*rx(n)*dt r(n) = -ri aa(n) = -ra*cv(nm1) bb(n) = ch(n) -aa(n) +rb*ck(n)/dzl cc(n) = 0.0d0 c*------- Jacobian enhancements c aa(n) = aa(n) -rb*0.25d0*(dk(nm1)+dk(n)) c bb(n) = bb(n) +rb*(dk(n) -0.25d0*(dk(nm1)+dk(n))) c*------- test well-posedness of tridiagonal solution c do i = 1, n c if (dabs(bb(i)) .eq. 0.d0) noluck = 2 c if (dabs(bb(i)) .lt. (dabs(aa(i)) + dabs(cc(i))) ) c + noluck = 2 c if (noluck .eq. 2) then c print *, 'poorly posed tridiag in water at i =', i c return c end if c end do c*------- ensure diagonal dominance of tridiagonal solution c do i = 1, n c dj = dabs(aa(i)) + dabs(cc(i)) - dabs(bb(i)) c if (dj .gt. 0.d0) then c bb(i) = bb(i) + dsign(dj, bb(i)) c end if c end do c*-------- call thomas algorithm to solve for soil water pressure call dgtsl (n, aa, bb, cc, r, info) if (info .gt. 0) then noluck = 2 return end if c*------- advance p(i) by Newton correction do i = 1, n c p(i) = p(i) + r(i) p(i) = p(i) + ur*r(i) c if(p(i) .gt. pmax) p(i) = pmax ht(i) = p(i) - z(i) end do if (zrj .le. 0.d0) then zr = ht(1) - 0.5d0*qe*dz/cv(0) else if (qinp .eq. zrj/dt-qe) then zr = ht(1) +0.5d0*dz*qinp/cv(0) else zr = ( zrj +ap*(zrp-zrj) -a2*qe*dt +2.d0*rb*cv(0)*ht(1) ) + /(1.d0+2.d0*rb*cv(0)) end if end if ht(0) = zr p(0) = zr c* --------- update soil water properties call wprop c*------- calculate et, root extraction and surface evaporation call etrx c*------- calculate residual mass balance (cm^3/cm^2) c if (zrp .le. 0.d0) then if (zrj .le. 0.d0) then qin = -qe zr = ht(1) - 0.5d0*qe*dz/cv(0) else c zr = ( zrj +ap*(zrp-zrj) -a2*qe*dt +2.d0*rb*cv(0)*ht(1) ) c + /(1.d0+2.d0*rb*cv(0)) c if (zr .lt. 0.d0) then c zr = 0.d0 c qin = zrp/dt -qe c else c qin = -2.d0*cv(0)*(ht(1)-zr)/dz c end if if (qinp .eq. zrj/dt-qe) then zr = ht(1) +0.5d0*dz*qinp/cv(0) qin = qinp else zr = ( zrj +ap*(zrp-zrj) -a2*qe*dt +2.d0*rb*cv(0)*ht(1) ) + /(1.d0+2.d0*rb*cv(0)) qin = -2.d0*cv(0)*(ht(1)-zr)/dz end if end if ht(0) = zr p(0) = zr ql = -ck(n)*(ht(np1)-ht(n))/dzl qn = qin - ql sdth = 0.d0 sresp = 0.d0 do i = 1, n c sdth = sdth + th(i)-thj(i) sdth = sdth + th(i)-thj(i) -ap*(thp(i)-thj(i)) qn = qn - rx(i)*dz sresp = sresp + dabs(r(i)*1d-2) end do c res1 = sdth*dz - (a1*qnp + a2*qn)*dt res1 = sdth*dz - a2*qn*dt res2 = sresp res = dsqrt(dabs(res1)*res2) c res = res1 if (inewt .eq. 1) then resl = res url = ur else if (dabs(res) .le. dabs(resl)) then ur = ur + 0.1d0*(urmax-ur) do i = 0, n hti(i) = ht(i) end do else if (ur .gt. 0.1d0) then ur = dmax1(0.1d0, 0.5d0*ur) do i = 0, n ht(i) = hti(i) p(i) = ht(i) + z(i) end do call wprop call etrx go to 10 end if end if write (*, '(1h+, i7, a, i7, 2f15.10, f6.3)') + inewt, ' NEWT IFC', ifc, res, dt, ur c if (dt .ge. 0.1d0) then if (dabs(res) .lt. 1.d-6) go to 70 c else c if (dabs(res) .lt. 1.d-5) go to 70 c end if resl = res 60 continue 70 continue if (dt .lt. 1.d-100) then noluck = 4 print *, 'dt goes to zero in subroutine water' return end if if (inewt .gt. maxit) then dt = 0.619*dt noluck = -1 if (dt .lt. 1.d-100) then print *, 'endless loop in subroutine water' noluck = 1 return end if return end if c*------- update cqin if (zrj .gt. 0.d0) then c cqin = cqin +dt*(a1*qinp +a2*qin) cqin = cqin +dt*a1*qinp end if if (zrp .gt. 0.d0) then cqin = cqin +dt*a2*qin end if write (*, '(1h+, I7, A, I7, 2f15.10, f6.3, f9.4)') + inewt, ' NEWT IFC', ifc, res, dt, ur, cqin c*------- update zrj, zrp if (rday .and. (zrj .gt. 0.d0) .and. (zr .le. 0.d0) + .and. (inewt .lt. maxit)) then dtfl = 1 end if zrp = zr zrj = zr c*------- calculate pore water velocity and direction arrays do i = 1, n if (i .lt. n) then dh = (ht(i)-ht(i+1))/dz if (dh .ge. 0.d0) then vd(i) = .true. v(i) = cv(i)*dh else vd(i) = .false. v(i) = cv(i)*dh end if else dh = (ht(n)-ht(np1))/(dzl+0.5d0*dz) if (dh .ge. 0.d0) then vd(n) = .true. else vd(n) = .false. end if v(n) = ck(n)*dh end if end do c* --------- calculate water extracted by roots from the root zone sum = 0.d0 do i = 1, n sum = sum + rx(i) end do plnw = plnw + (a1*srxp+a2*sum)*dt*dz cev = cev + qe*dt return end c*********************************************************************** c* this subroutine reads in weather data from a dataset in an input c* data file it calls the readit subroutine. it also calculates some of c* the constants required by hargreaves formula for estimating evapotr- c* anspiration, and calls the evapo subroutine to calculate the daily c* evapotranspiration. c*********************************************************************** subroutine weather(allok) implicit real*8(a-h,o-z) include 'come4.for' logical allok xlat = 36.d0 c* --------- call readit subroutine to read in weather c* ----------------- data from the dataset call readit(allok) if (allok) then c*------- roll over time series arrays at first of year if ((kyer(klm) .gt. kyer(klm-1)) .and. (klm .ne. 1)) then accpln=0.d0 c accevp=0.d0 do ib = -4,0,1 crdaly(ib) = crdaly(klm+ib-1) end do c csno3(1) = csno3(klm) et(1) = et(klm) kdoy(1) = kdoy(klm) kyer(1) = kyer(klm) pet(1) = pet(klm) runof(1) = runof(klm) tempmx(1) = tempmx(klm) tempmn(1) = tempmn(klm) xinfl(1) = xinfl(klm) klm = 1 lqtr = 0 end if c* ---------- set up the constants that are required for the c* -------------- calculation of the evapotranspiration rm = 3206.18d0 am = -27405.5d0 bm = 54.1896d0 cm = -0.045137d0 dm = 0.215321d-4 em = -0.462027d-8 fm = 2.41613d0 gm = 0.00121547d0 nday = kdoy(klm) c relative humidity and radiation on a horizontal surface for c eeq calculation in watbl2. (ashrea and handbook of meteorology c berry, bollay and beers 1945. p.293) tmaxk = 273.15d0 + tmax tmink = 273.15d0 + tmin tavek = 0.5d0*(tmaxk + tmink) tave = 0.5d0*(tmax + tmin) blt = arrhen (-6750.d0, tave) ps = (dexp((am + (bm*tmink) + (cm*tmink**2) + (dm*tmink**3) + + (em*tmink**4)) / ((fm*tmink) - (gm*tmink**2)))) * rm pv = (dexp((am + (bm *tavek)+(cm * tavek**2)+(dm * tavek**3)+ + (em * tavek**4)) / ((fm * tavek) - (gm * tavek**2)))) * rm rh = (ps / pv) * 100.d0 if (rh .lt. 100.d0) then xkt = 0.035d0 * ((100.d0- rh )**(.33d0)) else xkt = 0.d0 end if dec =(23.5d0*cos(((nday- + 173.d0)/365.25d0)*6.284d0))*3.142d0/180.0d0 hm = acos(-tan(xlat*0.01745d0) * tan(dec)) rad = (2880d0/3.142d0)*sin(xlat*0.01745d0) + *sin(dec)*(hm-tan(hm)) call evapo call jconvt c*------- set quarterly flag if end of quarter of year qtr = .false. if ((imm .eq. 3) .and. (idd .eq. 31)) qtr = .true. if ((imm .eq. 6) .and. (idd .eq. 30)) qtr = .true. if ((imm .eq. 9) .and. (idd .eq. 30)) qtr = .true. if ((imm .eq. 12) .and. (idd .eq. 31)) qtr = .true. endif return end c*************************************************************************** c* this subroutine provides the soil-water properties for each soil layer * c* in the soil profile, namely the hydraulic conductivity as a function of * c* water content or pressure and the water content-pressure relationship. * c* from the latter, the water capacity term is calculated. * c*************************************************************************** subroutine wprop implicit real*8(a-h,o-z) include 'come4.for' c*------- setup soil water properties call wpset(0, np1) c*------- Calculate Darcian-wieghted interblock conductivity means, cv do i = 0, nm1 zd = z(i+1) - z(i) cs = 0.5d0*(sc(i)+sc(i+1)) cb = 0.5d0*(bc(i)+bc(i+1)) c rht = 0.5d0*(thr(i)+thr(i+1)) c sht = 0.5d0*(ths(i)+ths(i+1)) ta = 0.5d0*(at(i)+at(i+1)) em = 0.5d0*(vm(i)+vm(i+1)) en = 0.5d0*(vn(i)+vn(i+1)) cv(i) = dkmh (ck(i), ck(i+1), p(i), p(i+1), cs, cb, + ta, em, en, zd, wv(i)) end do cv(n) = ck(n) return end double precision function dkmh (kk1, kk2, hh1, hh2, sc, bc, + at, vm, vn, dz, wv) implicit real*8 (a-h, k, o-z) implicit integer (i-j, l-n) c k1 = kk1/sc c k2 = kk2/sc k1 = kk1 k2 = kk2 h1 = hh1 h2 = hh2 tol = 1.d-6 if (dabs(k2-k1)/dsqrt(k1*k2) .le. tol) then dkmh = 0.5d0*(k1+k2) return end if if ((h1 .ge. 0.d0) .and. (h2 .ge. 0.d0)) then dkmh = 0.5d0*(k1+k2) return end if if (h1 .gt. -tol) h1 = -tol if (h2 .gt. -tol) h2 = -tol kl = dsqrt(k1*k2) if (kl .ge. sc) then dkmh = dink (-h1, -h2, sc, bc, at, vm, vn) return end if c go to 100 pl = at*((kl/sc)**(-vm/(vn*bc)) - 1.d0)**(1.d0/vm) pl1 = (pl/at)**vm u = bc*dz*pl1/(pl*(1.d0+pl1)) if (u .le. 0.0036d0) then wv = u/6.d0 else if (u .ge. 36.d0) then wv = (u-2.d0)/u else e = dexp(-u) wv = (u-2.d0+(u+2.d0)*e)/(u*(1.d0-e)) end if c 100 continue c wv = 0.1d0 kh = dink (-h1, -h2, sc, bc, at, vm, vn) c dkmh = 0.00944d0*(wv*k2 + (1.d0-wv)*kh) dkmh = (wv*k2 + (1.d0-wv)*kh) return end double precision function dink (p1, p2, sc, bc, at, vm, vn) implicit real*8 (a-h, k, o-z) dimension g(2), w(2) data (g(i), i=1,2) / 0.8611363116d0, 0.3399810436d0 / data (w(i), i=1,2) / 0.3478548451d0, 0.6521451549d0 / if (p1 .eq. p2) then dink = sc return end if dink = 0.d0 pm = .5d0*(p1+p2) pr = .5d0*(p2-p1) do i = 1, 2 pg = pr*g(i) dink = dink + w(i)* + ( ptc((pm+pg), sc, bc, at, vm, vn) + + ptc((pm-pg), sc, bc, at, vm, vn) ) end do dink = 0.5d0*dink return end double precision function ptc (ps, sc, bc, at, vm, vn) implicit real*8 (a-h, k, o-z) common / p2 / pd, be if (ps .le. 0.d0) then ptc = sc else ptc = sc/(1.d0+((ps/at)**vm))**(vn*bc/vm) end if return end c****************************************************************************** c this subroutine sets up soil water properties for the specified soil layer * c****************************************************************************** subroutine wpset(ll, lu) implicit real*8(a-h,o-z) include 'come4.for' ct = 0.00459d0/0.0022d0 - 1.d0 do 10 i=ll,lu pp = p(i) cmax = sc(i) tths = ths(i) tthr = thr(i) tbc = bc(i) tat = at(i) gn = vn(i) gm = vm(i) if (pp .ge. 0.0d0) then th(i) = tths ch(i) = 1.0d-20 c ch(i) = cmax*1.d-12 ck(i) = cmax dk(i) = 0.d0 else p1 = (-pp/tat)**gm p2 = p1/(pp*(1.d0+p1)) the = 1.0d0/(1.d0+p1)**(gn/gm) th(i) = tthr + the*(tths-tthr) ch(i) = -(tths-tthr)*the*gn*p2 ck(i) = cmax*the**tbc dk(i) = -tbc*gn*ck(i)*p2 if (ck(i) .le. 1.d-100) ck(i) = 1.d-100 if (ch(i) .le. 1.d-100) ch(i) = 1.d-100 c if (ch(i) .le. 1.d-100) ch(i) = ck(i)*1.d-12 end if tdif(i) = 0.00459d0/(1.d0+ct*dexp(-7.497d0*th(i))) + *3600.d0 10 continue return end c****************************************************************************** c* this function subprogram provides the rate coefficients for nitrification * c* as a function of soil water pressure (suction) for individual soil layers. * c****************************************************************************** c function zz1 (m, pp, ts) c****************************************************************************** c* this function subprogram provides the rate coefficients for denitrification* c* as a function of soil water content for individual soil layers. * c****************************************************************************** c function zz2 (m, wc, ts)