From 5e096283eb9ea892abed174aa5d28bbab600f896 Mon Sep 17 00:00:00 2001 From: crazyzlj Date: Wed, 22 Apr 2020 00:24:04 +0800 Subject: [PATCH] Update to 2012Rev591 --- VERSIONS | 2 +- src/allocate_parms.f | 8 +-- src/bmp_det_pond.f | 161 +++++++++++++++++++++--------------------- src/bmp_sand_filter.f | 10 +-- src/bmp_sed_pond.f | 12 ++-- src/bmpinit.f | 45 ++++++------ src/etact.f | 2 +- src/grow.f | 4 +- src/harvestop.f | 2 +- src/headout.f | 40 ----------- src/hrupond.f | 4 +- src/hydroinit.f | 7 +- src/irr_res.f | 2 +- src/main.f | 6 +- src/modparm.f | 1 + src/newtillmix.f | 2 +- src/nup.f | 4 +- src/plantmod.f | 2 +- src/pmeas.f | 4 +- src/pond.f | 2 +- src/print_hyd.f | 5 +- src/readbsn.f | 1 + src/readfile.f | 3 + src/readpnd.f | 1 + src/readrte.f | 16 ++--- src/readsol.f | 2 +- src/readsub.f | 15 ++-- src/res.f | 11 +-- src/rewind_init.f | 1 - src/routres.f | 14 ++-- src/sched_mgt.f | 21 +++--- src/simulate.f | 12 +++- src/soil_chem.f | 7 +- src/virtual.f | 22 +++--- src/watqual.f | 1 - src/watqual2.f | 60 ++++------------ src/writeaa.f | 2 +- 37 files changed, 231 insertions(+), 283 deletions(-) diff --git a/VERSIONS b/VERSIONS index fbfb2e1..2afa8af 100644 --- a/VERSIONS +++ b/VERSIONS @@ -1,3 +1,3 @@ VERSION_MAJOR 2012 -VERSION_MINOR 582 +VERSION_MINOR 591 VERSION_PATCH \ No newline at end of file diff --git a/src/allocate_parms.f b/src/allocate_parms.f index 8e0f98d..b4b4071 100644 --- a/src/allocate_parms.f +++ b/src/allocate_parms.f @@ -355,8 +355,8 @@ subroutine allocate_parms allocate (fcst_reg(msub)) allocate (harg_petco(msub)) ! allocate (hqd(nstep*3+1)) !! was 73, changed for urban - allocate (hqdsave(msub,nstep*2)) !! was 49, changed for urban -> changed to 2d array J.Jeong 4/17/2009 - allocate (hsdsave(msub,nstep*2)) !! J.Jeong 4/22/2009 + allocate (hqdsave(msub,nstep*4)) !! was 49, changed for urban -> changed to 2d array J.Jeong 4/17/2009 + allocate (hsdsave(msub,nstep*4)) !! J.Jeong 4/22/2009 !! allocate (hqd(73)) !! allocate (hqdsave(msub,49)) allocate (hru1(msub)) @@ -1184,7 +1184,7 @@ subroutine allocate_parms allocate (qdayout(mhru)) allocate (rch_dakm(mxsubch)) allocate (rchrg(mhru)) - allocate (rchrg_n(mhru)) + allocate (rchrg_n(mhru)) !! amount of nitrate getting to the shallow aquifer allocate (rchrg_dp(mhru)) allocate (re(mhru)) allocate (revapmn(mhru)) @@ -1509,7 +1509,7 @@ subroutine allocate_parms allocate (wshddayo(mstdo)) allocate (wshdmono(mstdo)) allocate (wshdyro(mstdo)) - allocate (fcstaao(nstep)) + allocate (fcstaao(16)) allocate (wpstaao(mpst,5)) allocate (wpstmono(mpst,5)) diff --git a/src/bmp_det_pond.f b/src/bmp_det_pond.f index 5a1257a..aa6a887 100644 --- a/src/bmp_det_pond.f +++ b/src/bmp_det_pond.f @@ -68,10 +68,14 @@ subroutine det_pond qpnd = dtp_ivol(sb) !m^3 sedpnd = dtp_ised(sb) !tons + ! Storage capacity under addon + qaddon = (dtp_addon(sb,1) / dtp_parm(sb)) ** 3. + & / (3.*ch_s(2,sb)) * pi !m3 + !! iterate for subdaily flow/sediment routing do ii=1,nstep - qout = 0.; qovmax = 0 + qout = 0.; qovmax = 0; depaddon = 0. qin = hhvaroute(2,ihout,ii) + qpnd !m^3 sedin = hhvaroute(3,ihout,ii) + sedpnd !tons if (qin>1e-6) then @@ -86,92 +90,87 @@ subroutine det_pond !! skip to next time step if no ponding occurs if (qdepth<=0.0001) cycle - !! Calculate weir outflow - do k=1,dtp_numstage(sb) + if (dtp_stagdis(sb)==0) then + !! Calculate weir outflow + do k = 1, dtp_numstage(sb) qstage = 0. - ! height to the top of addon from bottom - depaddon = depaddon + dtp_addon(sb,k) - if (k>1) depaddon = depaddon + dtp_depweir(sb,k-1) - ! volume below the current stage addon - qaddon = (depaddon / dtp_parm(sb)) ** 3. - & / (3.*ch_s(2,sb)) * pi !m3 - - !! Circular weir ? - if (dtp_weirtype(sb,k)==2) then - dtp_depweir(sb,k) = dtp_diaweir(sb,k) + dtp_addon(sb,k) - end if - - !! Fully submerged - if (qdepth>dtp_depweir(sb,k)) then - qdepth = qdepth - dtp_depweir(sb,k) - if (dtp_weirtype(sb,k)==1) then - watdepact = qdepth + dtp_depweir(sb,k) / 2 - warea = dtp_depweir(sb,k) * dtp_wrwid(sb,k) - else + + !! calculate weir discharge + + if (dtp_weirtype(sb,k)==2) then + !! Circular weir + dtp_depweir(sb,k) = dtp_diaweir(sb,k) + dtp_addon(sb,k) + + if (qdepth>dtp_depweir(sb,k)) then + !! Fully submerged + qdepth = qdepth - dtp_depweir(sb,k) watdepact = qdepth + dtp_diaweir(sb,k) / 2 warea = 3.14159 * dtp_diaweir(sb,k) ** 2 / 4. - end if - - !! Estimate discharge using orifice equation - qstage = dtp_cdis(sb,k) * 0.6 * warea * - & sqrt(19.6 * watdepact) !m3/s/unit - qstage = qstage * dtp_numweir(sb) * 60. * idt !m^3 - !Limit total outflow amount less than available water above addon - if(qin-qstageqaddon) qovmax = qin - qaddon - if(qstage>qovmax) qstage = qovmax - - !! Use stage-discharge relationship if available - if (dtp_stagdis(sb)==1) then - select case(dtp_reltype(sb)) - case(1) !! 1 is exponential function - qstage = dtp_coef1(sb) * exp(dtp_expont(sb) * qdepth) + - & dtp_intcept(sb) - case(2) !! 2 is Linear function - qstage = dtp_coef1(sb) * qdepth + dtp_intcept(sb) - case(3) !! 3 is logarthmic function - qstage = dtp_coef1(sb) * log(qdepth) + dtp_intcept(sb) - case(4) !! 4 is power function - qstage = dtp_coef1(sb) * (qdepth**3) + dtp_coef2(sb) * - & (qdepth**2) + dtp_coef3(sb) * qdepth + dtp_intcept(sb) - case(5) - qstage = dtp_coef1(sb)*(qdepth**dtp_expont(sb))+ - & dtp_intcept(sb) - end select - qstage = qstage * 60. * idt - end if !! end of stage-discharge calculation - qout = qout + qstage - exit - end if + qstage = qstage * dtp_numweir(sb) * 60. * idt !m^3 + + end if - - end do !! End of weir discharge calculations + qout = qout + qstage + end do + + !Limit total outflow amount less than available water above addon + if(qout>qin-qaddon) qout = max(qin - qaddon,0.) + + !! Flow over the emergency weir + watdepact = qdepth - (dtp_depweir(sb,1) + dtp_addon(sb,1)) + if (dtp_weirtype(sb,k)==1 .and. watdepact>0.) then + qstage = dtp_cdis(sb,k) * 1.84 * + & dtp_totwrwid(sb) * (watdepact ** 1.5) * 60. * idt !m3/s + qout = qout + qstage + end if - !! Check mss balance for flow + + else + !! Use stage-discharge relationship if available + if (dtp_stagdis(sb)==1) then + select case(dtp_reltype(sb)) + case(1) !! 1 is exponential function + qout = dtp_coef1(sb) * exp(dtp_expont(sb) * qdepth) + + & dtp_intcept(sb) + case(2) !! 2 is Linear function + qout = dtp_coef1(sb) * qdepth + dtp_intcept(sb) + case(3) !! 3 is logarthmic function + qout = dtp_coef1(sb) * log(qdepth) + dtp_intcept(sb) + case(4) !! 4 is power function + qout = dtp_coef1(sb) * (qdepth**3) + dtp_coef2(sb) * + & (qdepth**2) + dtp_coef3(sb) * qdepth + dtp_intcept(sb) + case(5) + qout = dtp_coef1(sb)*(qdepth**dtp_expont(sb))+ + & dtp_intcept(sb) + end select + qout = qout * 60. * idt + end if !! end of stage-discharge calculation + end if + + !! Check mass balance for flow if (qout>qin) then !no detention occurs qout = qin qpnd = 0. diff --git a/src/bmp_sand_filter.f b/src/bmp_sand_filter.f index 31772e6..f0baa4a 100644 --- a/src/bmp_sand_filter.f +++ b/src/bmp_sand_filter.f @@ -38,7 +38,7 @@ subroutine sand_filter(kk,flw,sed) & wetfsh,whd,sub_ha,dt,qcms,effct,effl,effg,effbr,vpipe,phead,hpnd, & tmpw,qloss,fsat,qpipe,mu,pipeflow,splw,hweir,tst,kb,qintns,qq, & qfiltr,sloss,spndconc,sedpnd,qpndi,qpnde,sedrmeff,sed_removed, - & sedconc,qevap + & sedconc,qevap,hrd real*8, dimension(:) :: qpnd(0:nstep),qsw(0:nstep),qin(0:nstep), & qout(0:nstep),fc(0:nstep),f(0:nstep) real, dimension(3,0:nstep), intent(inout) :: flw, sed @@ -164,8 +164,9 @@ subroutine sand_filter(kk,flw,sed) !soil water no more than saturation if (qsw(ii) > vfiltr) then - qout(ii) = ksat * dt / 1000. * tsa * ffsa !m3 - qsw(ii) = vfiltr + hrd = qsw(ii) / vfiltr + qout(ii) = ksat * hrd * dt / 1000. * tsa * ffsa !m3 + qsw(ii) = qsw(ii) - qout(ii) qpnd(ii) = qpndi - qsw(ii) - qout(ii) else if (qpnd(ii)>=qpnd(ii-1).and.qout(ii-1)<0.001) then @@ -243,7 +244,8 @@ subroutine sand_filter(kk,flw,sed) flw(3,ii) = qloss / (sub_ha *10000. - tsa) * 1000. !mm Endif - +! write(*,'(2i3,20f7.3)') iida, ii, qin(ii),qout(ii),qpnd(ii), +! & qsw(ii),qloss !-------------------------------------------------------------------------------------- ! TSS removal sloss = 0.; sedrmeff = 0. diff --git a/src/bmp_sed_pond.f b/src/bmp_sed_pond.f index bfbd9ee..8e4d17a 100644 --- a/src/bmp_sed_pond.f +++ b/src/bmp_sed_pond.f @@ -35,7 +35,7 @@ subroutine sed_pond(kk,flw,sed) real*8 :: tsa,mxvol,pdia,ksat,dp,sub_ha,mxh,hweir,phead,pipeflow real*8 :: qin,qout,qpnd,qpndi,sweir,spndconc,sedpnde,sedpndi,hpnd real*8 :: qweir, qtrns,qpipe,splw,sedconcweir,td,ksed,qevap - real, dimension(3,nstep), intent(inout) :: flw, sed + real, dimension(3,0:nstep), intent(inout) :: flw, sed sb = inum1 sub_ha = da_ha * sub_fr(sb) @@ -72,7 +72,7 @@ subroutine sed_pond(kk,flw,sed) If (hpnd > mxh) Then hweir = max(0.,hpnd - mxh) !water depth over weir crest, m !weir overflow, m^3 - qweir = 1.8 * (splw - 0.2 * hweir) * hweir ** 1.5 * + qweir = 3.33 * splw * hweir ** 1.5 * & idt * 60. hpnd = max(0.,(qpndi - qweir) / tsa) !m !overflow amount is no larger than surplus water above spillway height @@ -105,11 +105,11 @@ subroutine sed_pond(kk,flw,sed) !Outflow through orifice pipe hpnd = qpnd / tsa !m phead = max(0.,hpnd * 1000. - pdia / 2.) !mm - If (phead>pdia/2.) then +! If (phead>pdia/2.) then qpipe = pipeflow(pdia,phead) * idt *60. !m^3 - else - qpipe = qpnd * 0.1 - endif +! else +! qpipe = qout * 0.9 +! endif !update out flow, m^3 qout = qpipe diff --git a/src/bmpinit.f b/src/bmpinit.f index cf4ddae..93f14e8 100644 --- a/src/bmpinit.f +++ b/src/bmpinit.f @@ -75,9 +75,14 @@ subroutine bmpinit if (dtp_iyr(i)<=1000) dtp_iyr(i) = iyr if (dtp_evrsv(i)<=0) dtp_evrsv(i) = 0.1 if (dtp_numweir(i)<=0) dtp_numweir(i) = 1 - if (dtp_numstage(i)<=0) dtp_numstage(i) = 2 + if (dtp_numstage(i)<=0) dtp_numstage(i) = 1 + if (dtp_numstage(i)>1) then + do k=2,dtp_numstage(i) + if (dtp_weirtype(i,k)==1) dtp_addon(i,k) = 0. + end do + endif - !! Estimating emergency spillway volumes if not entered by user + !! Estimating design flow rate if not entered by user do k=1,dtp_numstage(i) if (dtp_flowrate(i,k)<=0.0) then dtp_flowrate(i,k) = 0.5 * 1000.0 * dtp_pcpret(i,k) @@ -92,10 +97,9 @@ subroutine bmpinit end if end do - !! Separate cumulative flow information to individual weir stages + !! Separate cumulative flow information to individual weir do k=2,dtp_numstage(i) - dtp_flowrate(i,k) = (dtp_flowrate(i,k) - & - sum(dtp_flowrate(i,1:k-1))) / dtp_numweir(i) + dtp_flowrate(i,k) = dtp_flowrate(i,k) / dtp_numweir(i) end do !!Estimate weir dimensions based on existing data @@ -108,12 +112,22 @@ subroutine bmpinit call est_weirdim(dtp_wdratio(i,k),dtp_flowrate(i,k) & ,dtp_wrwid(i,k),dtp_depweir(i,k),dtp_cdis(i,k)) end if - else - if (dtp_weirtype(i,k)==1) then !! reading user-entered data + else !! read user-entered data + if (dtp_weirtype(i,k)==1) then dtp_wrwid(i,k) = dtp_wdratio(i,k) * dtp_depweir(i,k) end if end if end do + + !! divide rectangular weirs into multiple single stage weirs + do k = 2, dtp_numstage(i) + dtp_addon(i,k) = dtp_addon(i,k-1) + dtp_depweir(i,k-1) + end do + + do k = dtp_numstage(i), 2, -1 + dtp_wrwid(i,k) = dtp_wrwid(i,k) - dtp_wrwid(i,k-1) + dtp_depweir(i,k-1) = dtp_depweir(i,k) + dtp_depweir(i,k-1) + end do end if !! Wet pond @@ -230,8 +244,8 @@ subroutine bmpinit wqv = hwq / 12. * sub_ha_urb(i) * sf_fr(i,k) * 107639.104167 !ft3 if (sf_dim(i,k)==0) then - write(77778,'(a46)') 'This SED-FIL size is automatically - & estimated' + write(77778,'(a46)') 'This SED-FIL size is automatically' + & // ' estimated based on WQV.' !Determine pond size automatically based on City of Austin's Design Guideline 1.6 if (sf_typ(i,k)==1.or.sf_typ(i,k)==3) then ! full scale or sedimentation basin only @@ -248,6 +262,7 @@ subroutine bmpinit ! wqv/(4.+1.33*4.) * 0.093 end if + sp_bpw(i,k) = 10. !m overflow weir width ft_pd(i,k) = 1524. !mm ft_dep(i,k) = 420. !mm ft_h(i,k) = 1200. !mm @@ -285,17 +300,7 @@ subroutine bmpinit !Orifice pipe for sand filter should be equal or larger than !sedimentation pond outlet pipe for full-type SedFils - if (ft_pd(i,k)1) sf_ptp(i,k) = 1 diff --git a/src/etact.f b/src/etact.f index 7b02f09..5de4a46 100644 --- a/src/etact.f +++ b/src/etact.f @@ -128,7 +128,7 @@ subroutine etact !! method is used to calculate surface runoff. The curve number methods !! take canopy effects into account in the equations. For either of the !! CN methods, canstor will always equal zero. - pet = pet - canstor(j) + pet = pet - canev if (pet < 0.) then canstor(j) = -pet canev = pet_day diff --git a/src/grow.f b/src/grow.f index 1a1e114..6ddd203 100644 --- a/src/grow.f +++ b/src/grow.f @@ -157,11 +157,11 @@ subroutine grow integer :: j real :: delg, par, ruedecl, beadj, reg, f, ff, deltalai - real :: laimax + real :: laimax, rto j = 0 j = ihru - + rto = 1. !! plant will not undergo stress if dormant if (idorm(j) == 1) return diff --git a/src/harvestop.f b/src/harvestop.f index ff8f16a..eb41035 100644 --- a/src/harvestop.f +++ b/src/harvestop.f @@ -418,7 +418,7 @@ subroutine harvestop if (laiday(j) < alai_min(idplt(j))) then !Sue laiday(j) = alai_min(idplt(j)) end if - if (phuacc(j) < .999) phuacc(j) = phuacc(j) * (1. - ff3) + phuacc(j) = phuacc(j) * (1. - ff3) rwt(j) = .4 - .2 * phuacc(j) else bio_ms(j) = 0. diff --git a/src/headout.f b/src/headout.f index e8cec7e..0e3535b 100644 --- a/src/headout.f +++ b/src/headout.f @@ -83,17 +83,6 @@ subroutine headout endif endif -!! write headings to VB interface HRU output file (hru.dat) - ilen = 10 - if (ipdvas(1) > 0) then - write (12,2000) (heds(ipdvas(j)), icols(ipdvas(j)), ilen, & - & j = 1, itots) - else - write (12,2000) (heds(j), icols(j), ilen, j = 1, mhruo) - endif - close (12) - - !! write headings to subbasin output file (output.sub) write (31,1000) prog, values(2), values(3), values(1), values(5), & & values(6), values(7) @@ -109,18 +98,6 @@ subroutine headout 1031 format (//6x,' SUB GIS MO DA YR AREAkm2',22(a10)) endif - -!! write headings to VB interface subbasin output file (sub.dat) - ilen = 10 - if (ipdvab(1) > 0) then - write (13,2000) (hedb(ipdvab(j)), icolb(ipdvab(j)), ilen, & - & j = 1, itotb) - else - write (13,2000) (hedb(j), icolb(j), ilen, j = 1, msubo) - endif - close (13) - - !! write headings to reach output file (output.rch) write (7,1000) prog, values(2), values(3), values(1), values(5), & & values(6), values(7) @@ -157,17 +134,6 @@ subroutine headout endif endif -!! write headings to VB interface reach output file (rch.dat) - ilen = 12 - if (ipdvar(1) > 0) then - write (11,2000) (hedr(ipdvar(j)), icolr(ipdvar(j)), ilen, & - & j = 1, itotr) - else - write (11,2000) (hedr(j), icolr(j), ilen, j = 1, mrcho) - endif - close (11) - - !! write headings to reservoir output file (output.rsv) write (8,1000) prog, values(2), values(3), values(1), values(5), & & values(6), values(7) @@ -180,12 +146,6 @@ subroutine headout write (22,1010) title write (22,1050) (hedrsv(j), j = 1, 41) end if - -!! write headings to VB interface reservoir output file (rsv.dat) - ilen = 12 - write (14,2000) (hedrsv(j), icolrsv(j-1), ilen, j = 2, 19) - close (14) - !! write headings to HRU impoundment output file (output.wtr) if (iwtr == 1) then diff --git a/src/hrupond.f b/src/hrupond.f index 93d5b1a..90ab546 100644 --- a/src/hrupond.f +++ b/src/hrupond.f @@ -113,7 +113,7 @@ subroutine hrupond !! calculate area of HRU covered by pond pndsa = 0. - pndsa = bp1(j) * pnd_vol(j) ** bp2(j) + pndsa = hru_fr(j) * bp1(j) * pnd_vol(j) ** bp2(j) !! calculate water flowing into pond for day pndflwi = qday + latq(j) @@ -179,7 +179,7 @@ subroutine hrupond call pond(j) !! compute water leaving pond - qday= qday + pndflwo / cnv +! qday = qday + pndflwo / cnv qdr(j) = qdr(j) + pndflwo / cnv !! compute sediment leaving pond diff --git a/src/hydroinit.f b/src/hydroinit.f index 33779ca..1dee670 100644 --- a/src/hydroinit.f +++ b/src/hydroinit.f @@ -203,13 +203,10 @@ subroutine hydroinit ql = q sumq = sumq + uh(isb,i) i = i + 1 - if (i>2*nstep) then - write(*,*) "Unit Hydrograph duration is too long" - exit - endif + if (i>3.*nstep) exit end do itb(isb) = i - 1 - do i = 1, itb(isb) + do i = 1, itb(isb) uh(isb,i) = uh(isb,i) / sumq end do endif diff --git a/src/irr_res.f b/src/irr_res.f index 3495bd6..e1af82a 100644 --- a/src/irr_res.f +++ b/src/irr_res.f @@ -106,7 +106,7 @@ subroutine irr_res !! check for timing of irrigation operation flag = 0 - flag = irra_flag(k) + flag = irr_flag(k) if (auto_wstr(k) > 0.) then if (wstrs_id(k) == 1 .and. strsw(k) < auto_wstr(k)) flag = 2 if (wstrs_id(k) == 2 .and. sol_sumfc(k) - sol_sw(k) > & diff --git a/src/main.f b/src/main.f index 3861d8a..f5b67b0 100644 --- a/src/main.f +++ b/src/main.f @@ -55,11 +55,11 @@ program main use parm implicit none - prog = "SWAT Nov 14 2012 VER 2012/Rev 582" + prog = "SWAT Apr 12 2013 VER 2012/Rev 591" write (*,1000) 1000 format(1x," SWAT2012 ",/, - & " Rev. 582 ",/, + & " Rev. 591 ",/, & " Soil & Water Assessment Tool ",/, & " PC Version ",/, & " Program reading from file.cio . . . executing",/) @@ -135,4 +135,4 @@ program main iscen=1 stop - end + end \ No newline at end of file diff --git a/src/modparm.f b/src/modparm.f index 7c80eb5..3fc740c 100644 --- a/src/modparm.f +++ b/src/modparm.f @@ -558,6 +558,7 @@ module parm real, dimension (:), allocatable :: det_san, det_sil, det_cla real, dimension (:), allocatable :: det_sag, det_lag real, dimension (:), allocatable :: tnylda, afrt_surface + real :: frt_surface real, dimension (:), allocatable :: auto_nyr, auto_napp real, dimension (:), allocatable :: manure_kg, auto_nstrs real, dimension (:,:), allocatable :: rcn_mo, rammo_mo diff --git a/src/newtillmix.f b/src/newtillmix.f index cca2321..c586708 100644 --- a/src/newtillmix.f +++ b/src/newtillmix.f @@ -7,7 +7,7 @@ subroutine newtillmix(jj,bmix) !! Mixing was extended to all layers !! A subroutine to simulate stimulation of organic matter decomposition was added !! March 2009: testing has been minimal and further adjustments are expected -!! use with caution and report anomalous results to akemanian@brc.tamus.edu and jeff.arnold@ars.usda.edu +!! use with caution !! ~ ~ ~ INCOMING VARIABLES ~ ~ ~ !! name |units |definition diff --git a/src/nup.f b/src/nup.f index 655842b..6799271 100644 --- a/src/nup.f +++ b/src/nup.f @@ -156,8 +156,8 @@ subroutine nup else xx = 1. end if - strsn = amax1(strsn, xx) - strsn = amin1(strsn, 1.) + strsn(j) = amax1(strsn(j), xx) + strsn(j) = amin1(strsn(j), 1.) end select return diff --git a/src/plantmod.f b/src/plantmod.f index b70cee5..c7effb6 100644 --- a/src/plantmod.f +++ b/src/plantmod.f @@ -60,7 +60,7 @@ subroutine plantmod sol_cov(j) = Max(sol_cov(j),0.) !! perform management operations - call operatn + ! call operatn !! compute plant water use and water stress !! compute actual plant transpiration diff --git a/src/pmeas.f b/src/pmeas.f index 717fc36..2fe8df1 100644 --- a/src/pmeas.f +++ b/src/pmeas.f @@ -92,6 +92,8 @@ subroutine pmeas ! allocate (hrmeas(mrg,25)) ! allocate (rhrbsb(24)) end if + + inum3sprev = 0 !! initialize variables for the day rmeas = 0. @@ -300,4 +302,4 @@ subroutine pmeas 5300 format (9x,a1) 5400 format (10x,"ERROR: Precipitation data dates do not match for", & & " simulation year: ",i4," and julian date: ",i3) - end + end \ No newline at end of file diff --git a/src/pond.f b/src/pond.f index dfe5bcc..49e80b8 100644 --- a/src/pond.f +++ b/src/pond.f @@ -154,7 +154,7 @@ subroutine pond(k) !! calculate water balance for day pndsa = 0. - pndsa = bp1(k) * pnd_vol(k) ** bp2(k) + pndsa = hru_fr(k) * bp1(k) * pnd_vol(k) ** bp2(k) pndev = 10. * evpnd(k) * pet_day * pndsa pndsep = pnd_k(k) * pndsa * 240. pndpcp = subp(k) * pndsa * 10. diff --git a/src/print_hyd.f b/src/print_hyd.f index 8b42ee4..a584d19 100644 --- a/src/print_hyd.f +++ b/src/print_hyd.f @@ -25,7 +25,7 @@ subroutine print_hyd use parm !! mauro/jerry whittaker hourly output file - if (iphr > 0) then + if (iphr > 0 .and. curyr > nyskip) then do ij = 1, nstep write (83,1000) iyr,i,ij,ihout,hhvaroute(2,ihout,ij) end do @@ -34,5 +34,4 @@ subroutine print_hyd !! end hourly codes return - end - + end \ No newline at end of file diff --git a/src/readbsn.f b/src/readbsn.f index 022f83a..7767f6d 100644 --- a/src/readbsn.f +++ b/src/readbsn.f @@ -570,6 +570,7 @@ subroutine readbsn !! set default values for undefined parameters ! if (ievent == 1) nstep = 24 + if (r2adj < 1.e-6) r2adj = 1. if (drain_co_bsn < 1.e-6) drain_co_bsn = 10. if (res_stlr_co < 1.e-6) res_stlr_co = .184 if (depimp_bsn < 1.e-6) depimp_bsn = 6000. diff --git a/src/readfile.f b/src/readfile.f index f9ec8cc..77f4d98 100644 --- a/src/readfile.f +++ b/src/readfile.f @@ -287,6 +287,9 @@ subroutine readfile read (101,5101) titldum read (101,*) iprint read (101,*) nyskip +!!! check that nyskip input is not greater or equal to nbyr input + if (nyskip >= nbyr) nyskip = nbyr - 1 + read (101,*) ilog read (101,*) iprp read (101,5101) titldum diff --git a/src/readpnd.f b/src/readpnd.f index 8199ec1..4ae883d 100644 --- a/src/readpnd.f +++ b/src/readpnd.f @@ -219,6 +219,7 @@ subroutine readpnd pndevcoeff = 0. wetevcoeff = 0. sub_ha = sub_km(i) * 100. + velsetlpnd = 0. do read (104,5100,iostat=eof) titldum diff --git a/src/readrte.f b/src/readrte.f index a51c36c..9f093ec 100644 --- a/src/readrte.f +++ b/src/readrte.f @@ -260,22 +260,22 @@ subroutine readrte !! An estimate of channel bank erodibility coefficient from jet test if it is not available !! Units of kd is (cm^3/N/s) !! Base on Hanson and Simon, 2001 - if (ch_bnk_kd(i) <= 1.e-6) then - if (tc_bnk(i) <= 1.e-6) then - ch_bnk_kd(i) = 0.2 + if (ch_bnk_kd(irch) <= 1.e-6) then + if (tc_bnk(irch) <= 1.e-6) then + ch_bnk_kd(irch) = 0.2 else - ch_bnk_kd(i) = 0.2 / sqrt(tc_bnk(i)) + ch_bnk_kd(irch) = 0.2 / sqrt(tc_bnk(irch)) end if end if !! An estimate of channel bed erodibility coefficient from jet test if it is not available !! Units of kd is (cm^3/N/s) !! Base on Hanson and Simon, 2001 - if (ch_bed_kd(i) <= 1.e-6) then - if (tc_bed(i) <= 1.e-6) then - ch_bed_kd(i) = 0.2 + if (ch_bed_kd(irch) <= 1.e-6) then + if (tc_bed(irch) <= 1.e-6) then + ch_bed_kd(irch) = 0.2 else - ch_bed_kd(i) = 0.2 / sqrt(tc_bed(i)) + ch_bed_kd(irch) = 0.2 / sqrt(tc_bed(irch)) end if end if diff --git a/src/readsol.f b/src/readsol.f index 980c385..971bbfe 100644 --- a/src/readsol.f +++ b/src/readsol.f @@ -121,7 +121,7 @@ subroutine readsol ! change below double subscripted sol_ec statement 1/27/09 when making septic changes read (107,5000,iostat=eof) (sol_ec(j,ihru), j = 1, nly) ! change below double subscripted sol_ec statement 1/27/09 when making septic changes - + !! MJW added rev 490 !!CaCo3 content (%) if (eof < 0) exit diff --git a/src/readsub.f b/src/readsub.f index 2f59fa5..fa3dcd0 100644 --- a/src/readsub.f +++ b/src/readsub.f @@ -378,13 +378,14 @@ subroutine readsub !! endif !! set default values -!! set default values - if (re(ihru) <= 0.) re(ihru) = re_bsn - if (sdrain(ihru) <= 0.) sdrain(ihru) = sdrain_bsn - if (drain_co(ihru) <= 0.) drain_co(ihru) = drain_co_bsn - if (pc(ihru) <= 0.) pc(ihru) = pc_bsn - if (latksatf(ihru) <= 0.) latksatf(ihru) = latksatf_bsn - if (sstmaxd(ihru) <= 0.) sstmaxd(ihru) = sstmaxd_bsn + do ihru = jj, hrutot(i) + if (re(ihru) <= 0.) re(ihru) = re_bsn + if (sdrain(ihru) <= 0.) sdrain(ihru) = sdrain_bsn + if (drain_co(ihru) <= 0.) drain_co(ihru) = drain_co_bsn + if (pc(ihru) <= 0.) pc(ihru) = pc_bsn + if (latksatf(ihru) <= 0.) latksatf(ihru) = latksatf_bsn + if (sstmaxd(ihru) <= 0.) sstmaxd(ihru) = sstmaxd_bsn + end do ! estimate drainage area for urban on-line bmps in square km !subdr_km(i) = subdr_km(i) + sub_km(i) diff --git a/src/res.f b/src/res.f index ab79e45..7137dca 100644 --- a/src/res.f +++ b/src/res.f @@ -155,11 +155,6 @@ subroutine res resflwo = res_out(jres,i_mo,curyr) !! This will override the measured outflow! This is just a check !! should really calibrate inflow or check res volumes - ndespill = ndtargr(nres) - if (ndespill <= 0.) ndespill = 10. - if (res_vol(jres) > res_evol(jres)) then - resflwo = resflwo+(res_vol(jres)-res_evol(jres))/ndespill - endif case (2) !! controlled outflow-target release targ = 0. @@ -208,6 +203,12 @@ subroutine res if (resflwo < oflowmn_fps(jres)) resflwo = oflowmn_fps(jres) end select + + ndespill = ndtargr(nres) + if (ndespill <= 0.) ndespill = 10. + if (res_vol(jres) > res_evol(jres)) then + resflwo = resflwo+(res_vol(jres)-res_evol(jres))/ndespill + endif !! if reservoir volume is zero if (res_vol(jres) < 0.001) then diff --git a/src/rewind_init.f b/src/rewind_init.f index 9f6b190..532c45e 100644 --- a/src/rewind_init.f +++ b/src/rewind_init.f @@ -139,7 +139,6 @@ subroutine rewind_init igen = 0 igen = ogen - !! HRU variables snoeb = 0. snoeb = orig_snoeb diff --git a/src/routres.f b/src/routres.f index fa08c5c..32fa69b 100644 --- a/src/routres.f +++ b/src/routres.f @@ -90,21 +90,21 @@ subroutine routres !! |reactions during month !! resoutm(8,:) |mg pst |pesticide lost from reservoir through !! |volatilization during month -!! resoutm(9,:) |mg pst |pesticide moving from water to sediment +!! resoutm(9,:) |mg pst |pesticide moving from water to sediment < !! |through settling during month -!! resoutm(10,:)|mg pst |pesticide moving from sediment to water +!! resoutm(10,:)|mg pst |pesticide moving from sediment to water < !! |through resuspension during month -!! resoutm(11,:)|mg pst |pesticide moving from water to sediment +!! resoutm(11,:)|mg pst |pesticide moving from water to sediment < !! |through diffusion during month -!! resoutm(12,:)|mg pst |pesticide lost from reservoir sediment layer +!! resoutm(12,:)|mg pst |pesticide lost from reservoir sediment layer < !! |through reactions during month -!! resoutm(13,:)|mg pst |pesticide lost from reservoir sediment layer +!! resoutm(13,:)|mg pst |pesticide lost from reservoir sediment layer < !! |through burial during month -!! resoutm(14,:)|mg pst |pesticide transported out of reservoir during +!! resoutm(14,:)|mg pst |pesticide transported out of reservoir during < !! |month !! resoutm(15,:)|mg pst/m^3 |pesticide concentration in reservoir water !! |during month -!! resoutm(16,:)|mg pst/m^3 |pesticide concentration in reservoir sediment +!! resoutm(16,:)|mg pst/m^3 |pesticide concentration in reservoir sediment < !! |layer during month !! resoutm(17,:)|m^3 H2O |evaporation from reservoir during month !! resoutm(18,:)|m^3 H2O |seepage from reservoir during month diff --git a/src/sched_mgt.f b/src/sched_mgt.f index 72e4683..c751e51 100644 --- a/src/sched_mgt.f +++ b/src/sched_mgt.f @@ -71,6 +71,8 @@ subroutine sched_mgt irr_flag(ihru) = 1 if (irrefm(ihru) < 1.e-6) irrefm(ihru)=1.0 + if (irr_sc(j) <= 0) irr_sc(j) = irrsc(j) + if (irr_no(j) <= 0) irr_no(j) = irrno(j) if (irr_no(j) <= 0) irr_no(j) = hru_sub(j) if (irr_sc(ihru) > 2) then !! reach and res flag ?? call irrsub @@ -81,7 +83,7 @@ subroutine sched_mgt * iida, " ", * "IRRIGATE", phubase(j), phuacc(j), sol_sw(j),bio_ms(j), * sol_rsd(1,j), sol_sumno3(j),sol_sumsolp(j),irramt(j), - * irrsc(j), irrno(j) + * irr_sc(j), irr_no(j) 1002 format (a5,1x,a4,3i6,2a15,7f10.2,10x,f10.2,70x,2i7) end if @@ -218,6 +220,7 @@ subroutine sched_mgt irr_asq(j) = mgt7op(nop(j),j) irr_sca(j) = mgt2iop(nop(j),j) irr_noa(j) = mgt10iop(nop(j),j) + if (irr_noa(j) <= 0) irr_noa(j) = irrno(j) if (irr_noa(j) <= 0) irr_noa(j) = hru_sub(j) if (wstrs_id(j) <= 0.) wstrs_id(j) = 1. if (irr_eff(j) > 1.) irr_eff(j) = 0. @@ -235,16 +238,18 @@ subroutine sched_mgt auto_nyr(j) = mgt6op(nop(j),j) if (auto_nyr(j) < 1.e-6) auto_nyr(j) = 350. auto_eff(j) = mgt7op(nop(j),j) - if (auto_eff(ihru) <= 0.) auto_eff(ihru) = 1.3 + if (auto_eff(j) <= 0.) auto_eff(j) = 1.3 afrt_surface(j) = mgt8op(nop(j),j) - if (afrt_surface(ihru) <= 1.e-6) afrt_surface(ihru) = .8 + if (afrt_surface(j) <= 1.e-6) afrt_surface(j) = .8 !! calculate tnylda for autofertilization ncrp = idplt(j) - if (hvsti(ncrp) < 1.) then - tnylda(j) = 350. * cnyld(ncrp) * bio_e(ncrp) - else - tnylda(j) = 1000. * cnyld(ncrp) * bio_e(ncrp) - endif + if (tnylda(j) < 1.e-6)tnylda(j)=150.*cnyld(ncrp)*bio_e(ncrp) + ! if (tnylda(j) < 1.e-6)tnylda(j)=350.*cnyld(ncrp)*bio_e(ncrp) + ! tnylda(j) = 350. * cnyld(ncrp) * bio_e(ncrp) + ! tnylda(j) = 350. * cnyld(ncrp) * bio_e(ncrp) + ! else + ! tnylda(j) = 1000. * cnyld(ncrp) * bio_e(ncrp) + ! endif case (12) !! street sweeping (only if iurban=2) diff --git a/src/simulate.f b/src/simulate.f index b62de73..09cc4ff 100644 --- a/src/simulate.f +++ b/src/simulate.f @@ -238,14 +238,14 @@ subroutine simulate call operatn dorm_flag = 0 endif - nop(ihru) = nop(ihru) + 1 + nop(ihru) = nop(ihru) + 1 if (nop(ihru) > nopmx(ihru)) then nop(ihru) = 1 end if phubase(ihru) = 0. - yr_skip = 0 + yr_skip(ihru) = 0 endif endif @@ -308,7 +308,13 @@ subroutine simulate end if phubase(j) = 0. - yr_skip = 0 + yr_skip(j) = 0 + endif + if (mgtop(nop(j),j) == 17) then + nop(j) = nop(j) + 1 + if (mgtop(nop(j),j) == 17) then + yr_skip(j) = 1 + end if endif end do diff --git a/src/soil_chem.f b/src/soil_chem.f index 988d6f8..f6d7bce 100644 --- a/src/soil_chem.f +++ b/src/soil_chem.f @@ -255,11 +255,10 @@ subroutine soil_chem !! Limit PSP range if (psp <.05) then psp = 0.05 - end if - if (psp > 0.9) then + else if (psp > 0.9) then psp = 0.9 end if - end if + end if sol_actp(j,i) = sol_solp(j,i) * (1. - psp) / psp @@ -401,7 +400,7 @@ subroutine soil_chem sol_LMN(j,i)=.1*sol_LMC(j,i) sol_LSC(j,i)=.42*sol_LS(j,i) sol_LSLC(j,i)=.8*sol_LSC(j,i) - sol_LSLNC(j,i)=.2*sol_LSC(j,i) + sol_LSLNC(j,i)=.2*sol_LSC(j,i) sol_LSN(j,i)=sol_LSC(j,i)/150. !sol_WOC(j,ihru)=sol_WOC(j,ihru)+sol_LSC(j,ihru)+sol_WLMC(j,ihru) sol_WOC(j,i)=sol_WOC(j,i)+sol_LSC(j,i)+sol_LMC(j,i) diff --git a/src/virtual.f b/src/virtual.f index 68e1993..f86a7ad 100644 --- a/src/virtual.f +++ b/src/virtual.f @@ -265,7 +265,7 @@ subroutine virtual integer :: j, sb, kk, ii real :: cnv, sub_ha, wtmp, baseflw, bf_fr,hr - real :: sub_hwyld(nstep), hqd(3*nstep), hsd(3*nstep),hqdtst(nstep) ! hqd, hsd locally defined. J.Jeong 4/26/2009 + real :: sub_hwyld(nstep), hqd(4*nstep), hsd(4*nstep),hqdtst(nstep) ! hqd, hsd locally defined. J.Jeong 4/26/2009 j = ihru sb = inum1 @@ -443,14 +443,14 @@ subroutine virtual if (ievent >= 2) then ! ! subdaily surface runoff, upland sediment for the subbasin - sub_ubnrunoff(sb,1:nstep) = sub_ubnrunoff(sb,1:nstep) - & / sub_fr(sb) - sub_ubntss(sb,1:nstep) = sub_ubntss(sb,1:nstep) / sub_fr(sb) - - sub_hhqd(sb,1:nstep) = sub_hhqd(sb,1:nstep) / sub_fr(sb) - sub_hhsedy(sb,1:nstep) = sub_hhsedy(sb,1:nstep) - & / sub_fr(sb) - sub_atmp(sb,1:nstep) = sub_atmp(sb,1:nstep) / sub_fr(sb) +! sub_ubnrunoff(sb,1:nstep) = sub_ubnrunoff(sb,1:nstep) +! & / sub_fr(sb) +! sub_ubntss(sb,1:nstep) = sub_ubntss(sb,1:nstep) / sub_fr(sb) +! +! sub_hhqd(sb,1:nstep) = sub_hhqd(sb,1:nstep) / sub_fr(sb) +! sub_hhsedy(sb,1:nstep) = sub_hhsedy(sb,1:nstep) +! & / sub_fr(sb) +! sub_atmp(sb,1:nstep) = sub_atmp(sb,1:nstep) / sub_fr(sb) !---------------------------------------------------- ! Simulate distributed urban BMPs in the subbasin @@ -561,7 +561,7 @@ subroutine virtual varoute(27,ihout) = sub_dlag(sb) !! detached lrg ag varoute(29,ihout) = sub_qd(sb) * sub_ha * 10. !! surface runoff varoute(30,ihout) = sub_latq(sb) * sub_ha * 10. !! lateral flow - varoute(30,ihout) = sub_tileq(sb) * sub_ha * 10. !! tile flow + varoute(31,ihout) = sub_tileq(sb) * sub_ha * 10. !! tile flow varoute(32,ihout) = sub_gwq(sb) * sub_ha * 10. !! groundwater flow !! varoute array has space for 33 different routing components @@ -632,8 +632,6 @@ subroutine virtual hhvaroute(20,ihout,ii) = 0. !!cmetal#1 hhvaroute(21,ihout,ii) = 0. !!cmetal#2 hhvaroute(22,ihout,ii) = 0. !!cmetal#3 - hhvaroute(23,ihout,ii) = baseflw - end if end do end if diff --git a/src/watqual.f b/src/watqual.f index 71ec2eb..157a0fd 100644 --- a/src/watqual.f +++ b/src/watqual.f @@ -431,7 +431,6 @@ subroutine watqual rch_cbod(jrch) = cbodcon - (yy + zz) * tday !!deoxygenation rate - rk1(jrch) = .5 coef = exp(-Theta(rk1(jrch),thrk1,wtmp) * tday) cbodrch = coef * cbodcon !!cbod rate loss due to settling diff --git a/src/watqual2.f b/src/watqual2.f index 8fe5ea7..562b383 100644 --- a/src/watqual2.f +++ b/src/watqual2.f @@ -246,11 +246,6 @@ subroutine watqual2 if (rtwtr / 86400.> 0.01.and.wtrin>0.001) then !! concentrations - - - - - !! initialize concentration of nutrient in reach wtrtot = 0. algcon = 0. @@ -263,10 +258,7 @@ subroutine watqual2 cbodcon = 0. o2con = 0. wtrtot = wtrin + rchwtr - - - - + algcon = algae(jrch) orgncon = organicn(jrch) nh3con = ammonian(jrch) @@ -278,8 +270,6 @@ subroutine watqual2 o2con = rch_dox(jrch) wtmp = wattemp(jrch) - - c write(104,*) 't',jrch,disoxin, wtrin, rch_dox(jrch) c o2con = (disoxin * wtrin + rch_dox(jrch) * rchwtr) / wtrtot @@ -330,9 +320,7 @@ subroutine watqual2 bc2mod = bc2(jrch) * cordo !! end O2 impact calculations - C tday is the calculation time step = 1 day - tday = 1.0 !! algal growth @@ -387,8 +375,7 @@ subroutine watqual2 gra = 0. endif end select - - + !! calculate algal biomass concentration at end of day !! (phytoplanktonic algae) !! QUAL2E equation III-2 @@ -418,9 +405,7 @@ subroutine watqual2 !! calculate dissolved oxygen concentration if reach at !! end of day QUAL2E section 3.6 equation III-28 - - - + uu = 0. vv = 0. ww = 0. @@ -452,9 +437,7 @@ subroutine watqual2 end if ddisox = o2con + (uu + vv - ww - xx - yy - zz) * tday o2proc=o2con-ddisox - if (ddisox < 0.00001) ddisox = 0.00001 - - + if (ddisox < 0.00001) ddisox = 0.00001 !! end oxygen calculations !! nitrogen calculations @@ -470,13 +453,12 @@ subroutine watqual2 dorgn = orgncon + (xx - yy - zz) * tday if (dorgn < 0.00001) dorgn = 0.00001 - !! calculate fraction of algal nitrogen uptake from ammonia !! pool QUAL2E equation III-18 f1 = 0. f1 = p_n * nh3con / (p_n * nh3con + (1. - p_n) * no3con + & & 1.e-6) - + !! calculate ammonia nitrogen concentration at end of day !! QUAL2E section 3.3.2 equation III-17 ww = 0. @@ -501,8 +483,6 @@ subroutine watqual2 dno2 = no2con + (yy - zz) * tday if (dno2 < 1.e-6) dno2 = 0. - - !! calculate nitrate concentration at end of day !! QUAL2E section 3.3.4 equation III-20 yy = 0. @@ -538,17 +518,8 @@ subroutine watqual2 dsolp = 0. dsolp = solpcon + (xx + yy - zz) * tday if (dsolp < 1.e-6) dsolp = 0. - - - - !! end phosphorus calculations - - - - - wtrtot = wtrin + rchwtr !! initialize inflow concentrations chlin = 0. @@ -577,7 +548,6 @@ subroutine watqual2 heatin = varoute(1,inum2) * (1. - rnum1) end if - wattemp(jrch) =(heatin * wtrin + wtmp * rchwtr) / wtrtot algae(jrch) = (algin * wtrin + dalgae * rchwtr) / wtrtot @@ -617,16 +587,16 @@ subroutine watqual2 disolvp(jrch) = 0.0 rch_cbod(jrch) = 0.0 rch_dox(jrch) = 0.0 - dalgae = 0.0 - dchla = 0.0 - dorgn = 0.0 - dnh4 = 0.0 - dno2 = 0.0 - dno3 = 0.0 - dorgp= 0.0 - dsolp = 0.0 - dbod = 0.0 - ddisox = 0.0 + dalgae = 0.0 + dchla = 0.0 + dorgn = 0.0 + dnh4 = 0.0 + dno2 = 0.0 + dno3 = 0.0 + dorgp= 0.0 + dsolp = 0.0 + dbod = 0.0 + ddisox = 0.0 soxy = 0.0 endif diff --git a/src/writeaa.f b/src/writeaa.f index cf9a1ef..0d4bd90 100644 --- a/src/writeaa.f +++ b/src/writeaa.f @@ -349,7 +349,7 @@ subroutine writeaa aairr = aairr / yrs do j = 1, nhru do nicr = 1, mcrhru(j) - if (ncrops(nicr,j) > 0) then + if (ncrops(nicr,j) > 0) then yldn(nicr,j) = yldkg(nicr,j) / ncrops(nicr,j) bio_aahv(nicr,j) = bio_hv(nicr,j) / ncrops(nicr,j) else