diff --git a/VERSIONS b/VERSIONS index e7ef3a6..0f0c14a 100644 --- a/VERSIONS +++ b/VERSIONS @@ -1,3 +1,3 @@ VERSION_MAJOR 2012 -VERSION_MINOR 670 +VERSION_MINOR 681 VERSION_PATCH \ No newline at end of file diff --git a/src/apex_day.f b/src/apex_day.f index 45de929..0d97658 100644 --- a/src/apex_day.f +++ b/src/apex_day.f @@ -119,7 +119,7 @@ subroutine apex_day !! changes for Mike Winchell if apex.swt files has a different start date than the SWAT run if (ifirsta(inum1) == 1) then do - read (112+inum1,*) idapa(inum1), iypa(inum1), flodaya(inum1), + read (400+inum1,*) idapa(inum1), iypa(inum1), flodaya(inum1), & seddaya(inum1), orgndaya(inum1), orgpdaya(inum1), no3daya(inum1), & minpdaya(inum1) if(idapa(inum1) == id1 .and. iypa(inum1) == iyr) exit @@ -149,7 +149,7 @@ subroutine apex_day varoute(21,ihout) = 0.0 varoute(22,ihout) = 0.0 if (curyr /= nbyr .and. iida /= idal) then - read (112+inum1,*) idapa(inum1), iypa(inum1), flodaya(inum1), + read (400+inum1,*) idapa(inum1), iypa(inum1), flodaya(inum1), & seddaya(inum1), orgndaya(inum1), orgpdaya(inum1), no3daya(inum1), & minpdaya(inum1) endif diff --git a/src/carbon_new.f b/src/carbon_new.f index 700ce1c..d60055d 100644 --- a/src/carbon_new.f +++ b/src/carbon_new.f @@ -238,7 +238,7 @@ real*8 Function fnetmin(poold, R1, R2, hc, dummy, poolm, xinorg, cc1) else NPsoil = (sol_mass * sol_n(k,j) / 100.)/ sol_orgp(k,j) end if - sol_orgp(k,j) = dmin1(sol_orgp(k,j), 25.) + !sol_orgp(k,j) = dmin1(sol_orgp(k,j), 25.) !commented 110118 Armen K CPsoil = CNsoil * NPsoil if (sol_rsd(k,j) > 0.00001) then diff --git a/src/command.f b/src/command.f index a661206..846234c 100644 --- a/src/command.f +++ b/src/command.f @@ -169,7 +169,7 @@ subroutine command case (13) call apex_day case (14) - call saveconc + if (curyr > nyskip) call saveconc case (17) call routeunit call sumhyd diff --git a/src/filter.f b/src/filter.f index 29d868f..b0cdee1 100644 --- a/src/filter.f +++ b/src/filter.f @@ -313,7 +313,7 @@ subroutine filter if (remove1 > 100.) remove1 = 100. if (remove1 < 0.) remove1 = 0. - remove21 = 29.3 + 0.51 * surq_remove2 + remove2 = 29.3 + 0.51 * surq_remove2 if (remove2 > 100.) remove2 = 100. if (remove2 < 0.) remove2 = 0. diff --git a/src/hhnoqual.f b/src/hhnoqual.f index 4564054..ce8b484 100644 --- a/src/hhnoqual.f +++ b/src/hhnoqual.f @@ -157,7 +157,8 @@ subroutine hhnoqual cbodcon = 0. o2con = 0. wtrtot = wtrin + hrchwtr(ii) - if (ii == 1) then + if (wtrtot > 0.01) then + if (ii == 1) then algcon = (algin * wtrin + algae(jrch) * hrchwtr(ii)) / wtrtot orgncon = (orgnin * wtrin + organicn(jrch) * hrchwtr(ii)) & / wtrtot @@ -185,9 +186,10 @@ subroutine hhnoqual solpcon = (dispin * wtrin + hsolp(ii-1) * hrchwtr(ii)) / wtrtot cbodcon = (cbodin * wtrin + hbod(ii-1) * hrchwtr(ii)) / wtrtot o2con = (disoxin * wtrin + hdisox(ii-1) * hrchwtr(ii)) / wtrtot + end if end if - if (algcon < 1.e-6) algon = 0.0 + if (algcon < 1.e-6) algcon = 0.0 if (orgncon < 1.e-6) orgncon = 0.0 if (nh3con < 1.e-6) nh3con = 0.0 if (no2con < 1.e-6) no2con = 0.0 diff --git a/src/hhwatqual.f b/src/hhwatqual.f index 9d8c56d..9a4a3f0 100644 --- a/src/hhwatqual.f +++ b/src/hhwatqual.f @@ -241,7 +241,7 @@ subroutine hhwatqual wtrin = 0. wtrin = hhvaroute(2,inum2,ii) * (1. - rnum1) - if (hrtwtr(ii) / (idt * 60.) > 0.01) then + if (hrtwtr(ii) / (idt * 60.) > 0.01.and.hdepth(ii) > 0.01) then !! concentrations !! initialize inflow concentrations chlin = 0. @@ -275,7 +275,7 @@ subroutine hhwatqual if (ammoin < 1.e-6) ammoin = 0.0 if (nitritin < 1.e-6) nitritin = 0.0 if (nitratin < 1.e-6) nitratin = 0.0 - if (orgnpin < 1.e-6) orgnpin = 0.0 + if (orgpin < 1.e-6) orgpin = 0.0 if (dispin < 1.e-6) dispin = 0.0 if (cbodin < 1.e-6) cbodin = 0.0 if (disoxin < 1.e-6) disoxin = 0.0 @@ -292,7 +292,8 @@ subroutine hhwatqual cbodcon = 0. o2con = 0. wtrtot = wtrin + hrchwtr(ii) - if (ii == 1) then + if (wtrtot > 0.01) then + if (ii == 1) then algcon = (algin * wtrin + algae(jrch) * hrchwtr(ii)) / wtrtot orgncon = (orgnin * wtrin + organicn(jrch) * hrchwtr(ii)) & / wtrtot @@ -320,8 +321,8 @@ subroutine hhwatqual solpcon = (dispin * wtrin + hsolp(ii-1) * hrchwtr(ii)) / wtrtot cbodcon = (cbodin * wtrin + hbod(ii-1) * hrchwtr(ii)) / wtrtot o2con = (disoxin * wtrin + hdisox(ii-1) * hrchwtr(ii)) / wtrtot - end if - + end if + endif if (algcon < 1.e-6) algcon = 0.0 if (orgncon < 1.e-6) orgncon = 0.0 if (nh3con < 1.e-6) nh3con = 0.0 diff --git a/src/hmeas.f b/src/hmeas.f index f6d1ac0..d038979 100644 --- a/src/hmeas.f +++ b/src/hmeas.f @@ -103,8 +103,9 @@ subroutine hmeas end do return -! 5200 format (7x,300f8.3) -! 5300 format (i4,i3,300f8.3) - 5200 format (7x,1800f8.3) - 5300 format (i4,i3,1800f8.3) + + !5200 format (7x,1800f8.3) + !5300 format (i4,i3,1800f8.3) + 5200 format (7x,1900f8.3) !! for Pouya + 5300 format (i4,i3,1900f8.3) !! for Pouya end \ No newline at end of file diff --git a/src/main.f b/src/main.f index 86bbf81..2cd238e 100644 --- a/src/main.f +++ b/src/main.f @@ -47,10 +47,10 @@ program main use parm implicit none - prog = "SWAT Sep 7 VER 2018/Rev 670" + prog = "SWAT May 26 VER 2020/Rev 681" write (*,1000) - 1000 format(1x," SWAT2018 ",/, - & " Rev. 670 ",/, + 1000 format(1x," SWAT2020 ",/, + & " Rev. 681 ",/, & " Soil & Water Assessment Tool ",/, & " PC Version ",/, & " Program reading from file.cio . . . executing",/) @@ -77,6 +77,8 @@ program main call std2 call openwth call headout + + !call sw_init !! convert integer to string for output.mgt file subnum = "" diff --git a/src/nrain.f b/src/nrain.f index 724ea60..7637827 100644 --- a/src/nrain.f +++ b/src/nrain.f @@ -69,7 +69,7 @@ subroutine nrain no3pcp = .01 * rcn_d(hru_sub(j)) * precipday sol_nh3(2,j) = sol_nh3(2,j) + nh3pcp + & drydep_nh4_d(hru_sub(j)) - sol_no3(2,j) = sol_no3(2,j) + no2pcp + + sol_no3(2,j) = sol_no3(2,j) + no3pcp + & drydep_no3_d(hru_sub(j)) end select diff --git a/src/openwth.f b/src/openwth.f index ff8cd65..9999982 100644 --- a/src/openwth.f +++ b/src/openwth.f @@ -93,19 +93,19 @@ subroutine openwth if (slrfile /= ' ') then !! open (137,file=slrfile,recl=800) - open (137,file=slrfile,recl=15000) + open (137,file=slrfile,recl=16000) read (137,5000) titldum end if if (rhfile /= ' ') then !! open (138,file=rhfile,recl=800) - open (138,file=rhfile,recl=15000) + open (138,file=rhfile,recl=16000) read (138,5000) titldum end if if (wndfile /= ' ') then !! open (139,file=wndfile,recl=800) - open (139,file=wndfile,recl=15000) + open (139,file=wndfile,recl=16000) read (139,5000) titldum end if diff --git a/src/potholehr.f b/src/potholehr.f index 22fc45c..4fd7d40 100644 --- a/src/potholehr.f +++ b/src/potholehr.f @@ -444,7 +444,7 @@ subroutine potholehr() if (iprint==3) then write (125,2000)i,j,curyr,k,pot_vol(j),potsa(j),spillo,potsep, & potev,sol_sw(j),potpcpmm,potflwi(j) / cnv, - & potsedi(j) / hru_ha(j),potflow,potsedo / hru_ha(j) + & potsedi(j) / hru_ha(j),potflwo,potsedo / hru_ha(j) endif if (pot_vol(j) > 1.e-6) then diff --git a/src/readatmodep.f b/src/readatmodep.f index f11aa82..681927d 100644 --- a/src/readatmodep.f +++ b/src/readatmodep.f @@ -43,7 +43,7 @@ subroutine readatmodep if (atmofile /= ' ') then open (127,file=atmofile) do iii = 1, 5 - read (127,*) titldum + read (127,5101) titldum end do else !! no filename present in file.cio - set defaults @@ -89,6 +89,6 @@ subroutine readatmodep !1001 format (3i6) !1002 format (1200f10.3) !1000 format (8x,4f10.3) -!5101 format (a80) +5101 format (a80) return end \ No newline at end of file diff --git a/src/readfig.f b/src/readfig.f index fcae6d2..1ac8de8 100644 --- a/src/readfig.f +++ b/src/readfig.f @@ -325,9 +325,9 @@ subroutine readfig call caps(apex_in) ! i = 0 ! i = inum1s(idum) - open (112+inum1s(idum),file=apex_in,recl=350) - do ii = 1, 10 - read (112+inum1s(idum),5200) titldum + open (400+inum1s(idum),file=apex_in,recl=350) + do ii = 1, 9 + read (400+inum1s(idum),5200) titldum end do !! code to read from apex output file diff --git a/src/readops.f b/src/readops.f index 10a1209..60c1673 100644 --- a/src/readops.f +++ b/src/readops.f @@ -95,7 +95,7 @@ subroutine readops do mon = 0 day = 0 - year = 0. + iyear = 0. mgt_op = 0 mgt1i = 0 mgt2i = 0 @@ -203,7 +203,7 @@ subroutine readops ro_bmp_flos(iops,ihru) = mgt10 !! Flow ro_bmp_seds(iops,ihru) = mgt11 !! Sediment ro_bmp_pps(iops,ihru) = mgt12 !! Particulate P - ro_bmp_sps(iops,ihru) = mg13 !! Soluble P + ro_bmp_sps(iops,ihru) = mgt13 !! Soluble P ro_bmp_pns(iops,ihru) = mgt14 !! Particulate N ro_bmp_sns(iops,ihru) = mgt15 !! Soluble N ro_bmp_bacs(iops,ihru) = mgt16 !! Bacteria diff --git a/src/readpnd.f b/src/readpnd.f index 3439cfd..2fcbf68 100644 --- a/src/readpnd.f +++ b/src/readpnd.f @@ -276,7 +276,7 @@ subroutine readpnd read (104,5100,iostat=eof) titldum if (eof < 0) exit if (titldum == ' '.or.titldum == 'Inputs used in')then - vselsetlpnd = 10.0 + velsetlpnd = 10.0 else backspace 104 read (104,*,iostat=eof) pnd_d50 diff --git a/src/rechour.f b/src/rechour.f index 909ee70..ba02942 100644 --- a/src/rechour.f +++ b/src/rechour.f @@ -195,6 +195,9 @@ subroutine rechour varoute(27,ihout) = varoute(27,ihout) + sedhr * 0. ! lag varoute(28,ihout) = 0. ! gravel + + QHY(ii,ihout,IHX(1)) = hhvaroute(2,ihout,ii) / (dthy * 3600.) !added by Jaehak for subdaily routing, 2019 + end do return diff --git a/src/route.f b/src/route.f index fd55355..5d382d4 100644 --- a/src/route.f +++ b/src/route.f @@ -136,6 +136,7 @@ subroutine route if (ievent > 0) then do ii = 1, nstep hrtwtr(ii) = hrtwtr(ii) + qdbank / dfloat(nstep) + hsdti(ii) = max(0.,hrtwtr(ii) / (dthy * 3600.)) !m3 -> m3/s end do end if diff --git a/src/rtday.f b/src/rtday.f index cf412e6..df16c14 100644 --- a/src/rtday.f +++ b/src/rtday.f @@ -159,13 +159,9 @@ subroutine rtday !! calculate volume of water leaving reach on day scoef = 0. rtwtr = 0. - !scoef = 2. * det / (2. * rttime + det) scoef = det / (rttime + det) if (scoef > 1.) scoef = 1. rtwtr = scoef * (wtrin + rchstor(jrch)) - !new storage coefficient replacement - rtwtr = vc * rcharea * 86400. - rtwtr = dmin1 (rtwtr, wtrin) !! calculate amount of water in channel at end of day rchstor(jrch) = rchstor(jrch) + wtrin - rtwtr diff --git a/src/rthmusk.f b/src/rthmusk.f index 8495fd9..be4f3a7 100644 --- a/src/rthmusk.f +++ b/src/rthmusk.f @@ -124,6 +124,7 @@ subroutine rthmusk jrch = 0 jrch = inum1 + hrtevp = 0 !! Compute storage time constant for reach xkm = 0. diff --git a/src/rthvsc.f90 b/src/rthvsc.f90 index da371fa..c9c5c61 100644 --- a/src/rthvsc.f90 +++ b/src/rthvsc.f90 @@ -66,7 +66,6 @@ subroutine rthvsc() !! inhyd |none |inflow hydrograph storage location number !! jrch |none |reach number !! p |m |wetted perimeter -!! scoef |none |storage coefficient !! nstep |none |No. of steps in a day (depends on model operational time step) !! topw |m |width of channel at water level !! vol |m^3 H2O |volume of water in reach @@ -81,22 +80,20 @@ subroutine rthvsc() use parm - integer :: jrch, ii, inhyd,j,l,cday,istep - real*8 :: wtrin, c, p, scoef - real*8 :: vol, topw,pcpday - real*8, dimension(nstep*5) :: QMS, QMSI,pcp + integer :: jrch, ii, inhyd,j,l + real :: wtrin, c, p + real :: vol, topw + real*8, dimension(nstep*100) :: QMS, QMSI real*8 :: ai, aii, ao, cbw, chw, fpw, g1, qi2, sss, xflo, zch, zi, zii, zo QMS = 0. QMSI = 0. hrtevp = 0. hrttlc = 0. - pcp(1:nstep) = precipdt(1:nstep) - pcpday=precipday jrch = inum1 inhyd = inum2 - + if (sum(QHY(:,inhyd,IHX(1)))==0.) return @@ -105,9 +102,10 @@ subroutine rthvsc() ii = ii - 1 DO K = 1, nstep+1 !number of time steps a day ii = ii + 1 - QMSI(ii) = QHY(K,inhyd,IHX(J)) !inflow from upstream during dt for the day, m^3 + QMSI(ii) = QHY(K,inhyd,IHX(J)) !inflow from upstream during dt for the day, m^3/s END DO - END DO + END DO + if(sum(QMSI(1:nstep))<=0.01) then if (sum(QHY(1:nstep,ihout,IHX(1)))<=0.01) then !no new inflow or outflow from previous days @@ -129,13 +127,8 @@ subroutine rthvsc() QI1 = QHY(1,inhyd,IHX(1)) ! m3/s !inflow at timestep 1 of the current day QO1 = QHY(1,ihout,IHX(1)) !outflow at timestep 1 of the current day that is predicted previous day - - !skip time steps with no inflow - DO L = 1, NHY(inhyd) - IF (QMSI(L)>0.) EXIT - QMS(L) = 0. - END DO - + !QMSI(1) = QMSI(1) + rchstor(jrch) / (dthy * 3600.) + ZCH = phi(7,jrch) !channel depth,m CBW = phi(6,jrch) ! channel bottom width, m XL3 = ch_l2(jrch) / 3.6 @@ -146,30 +139,32 @@ subroutine rthvsc() ADI = 0. STHY=0. - ! ROUTE UNTIL OUTFLOW IS GREATER THAN 0. - DO ii = L,NHY(inhyd) !routing starts when inflow > 0. - I1 = ii - 1 - QI2 = QMSI(ii) !inflow at the current timestep - ADI = ADI + QI2 - CALL HQDAV(AI,CBW,QI2,SSS,ZCH,ZI,CHW,FPW,jrch) !INPUT:CBW,QI2,SSS,ZCH output: AI,ZI,CHW,FPW - DD = SQRT((XLS+ZI)/XLS) - XFLO = QI2 / DD - CALL HQDAV(AII,CBW,XFLO,SSS,ZCH,ZII,CHW,FPW,jrch) - SMO = 2. * ADI - QI2 - XFLO = SMO - XLT * AII - !O2 = SUM OF INFLOW - STORAGE - IF(XFLO > 0.)EXIT - QMS(ii) = 0. - STHY = STHY + QI2 - ZOO = 0. - AOO = 0. - V = .5 * QMSI(ii) / AII !m/s - TT = XL3 / (V + 1.E-5) !hr - CVSC = MIN(.99,2.*dthy/(2.*TT+dthy)) - QI1 = QI2 - !WRITE(KW(1),6)NIT,TM,ZII,ZOO,AII,AOO,V,TT,CVSC,STHY,QI2,QO2 - END DO - ISM=ii + ! ROUTE UNTIL OUTFLOW IS GREATER THAN 0. + IF (rchstor(jrch)>0.001.and.QO1>0.0001) THEN + ii = 1 + QI2 = QI1 + QI1 = QHY(nstep,inhyd,IHX(4)) !Inflow during the previous time step (last time step of previous day), m3/s + STHY = rchstor(jrch) / (dthy * 3600.) + ELSE + + DO ii = 1, NHY(jrch) !routing starts when inflow > 0. + I1 = ii - 1 + QI2 = QMSI(ii) !inflow at the current timestep; m3/s + ADI = ADI + QI2 + CALL HQDAV(AI,CBW,QI2,SSS,ZCH,ZI,CHW,FPW,jrch) !INPUT:CBW,QI2,SSS,ZCH output: AI,ZI,CHW,FPW + DD = SQRT((XLS+ZI)/XLS) + XFLO = QI2 / DD + CALL HQDAV(AII,CBW,XFLO,SSS,ZCH,ZII,CHW,FPW,jrch) + SMO = 2. * ADI - QI2 + XFLO = SMO - XLT * AII + !O2 = SUM OF INFLOW - STORAGE + IF(XFLO > 0.)EXIT + QMS(ii) = 0. + STHY = STHY + QI2 + QI1 = QI2 + END DO + + END IF IIY=0 !Calculate discharge rate using variable storage coefficient @@ -187,31 +182,31 @@ subroutine rthvsc() IIY = 1 END IF DO IT = 1, 20 - CALL HQDAV(AO,CBW,G1,SSS,ZCH,ZO,CHW,FPW,jrch) - XX = MAX(.75,(XLS+ZI-ZO)/XLS) - DD = SQRT(XX) - V = MAX(.05,(QI2+G1)*DD/(AI+AO)) - T2 = XL3 / V - TT = .5 * (T1 + T2) - C = MIN(.99,2.*DTHY/(2.*TT+DTHY)) - Q1 = C * SIA !Q1 is outflow m3/s Eq(15) in Jeong et al. (2014) - GQ = Q1 - G1 - IF (abs(GQ/G1) < .001) EXIT - IF (GQ > 0.)THEN - GL = G1 - ELSE - GB = G1 - END IF - G1 = .5 * (GB + GL) - IF (GB-GL < .001)EXIT + CALL HQDAV(AO,CBW,G1,SSS,ZCH,ZO,CHW,FPW,jrch) + XX = MAX(.75,(XLS+ZI-ZO)/XLS) + DD = SQRT(XX) + V = MAX(.05,(QI2+G1)*DD/(AI+AO)) + T2 = XL3 / V + TT = .5 * (T1 + T2) + C = MIN(.99,2.*DTHY/(2.*TT+DTHY)) + Q1 = C * SIA !Q1 is outflow m3/s Eq(15) in Jeong et al. (2014) + GQ = Q1 - G1 + IF (abs(GQ/(G1+0.0001)) < .001) EXIT + IF (GQ > 0.)THEN + GL = G1 + ELSE + GB = G1 + END IF + G1 = .5 * (GB + GL)+0.0001 + IF (GB-GL < .001)EXIT END DO - QO2 = MAX(.01,G1) + QO2 = MAX(.0001,G1) QMS(ii) = QO2 STHY = SIA - QO2 T1 = T2 - IF(ii > NHY(inhyd)) THEN - QMSI(ii) = .9 * QMSI(I1) - IF(QMS(I1) <= 1.) EXIT + IF(ii > NHY(jrch)) THEN + QMSI(ii) = .9 * QMSI(I1) + IF(QMS(I1) <= 1.) EXIT END IF ii = ii + 1 QI1 = QI2 @@ -228,8 +223,8 @@ subroutine rthvsc() DO J = 2, 4 ii = ii - 1 DO K = 1, nstep+1 - ii = ii + 1 - QHY(K,ihout,IHX(J)) = QMS(ii) + ii = ii + 1 + QHY(K,ihout,IHX(J)) = QMS(ii) END DO END DO diff --git a/src/rtmusk.f b/src/rtmusk.f index a316976..17ed6ce 100644 --- a/src/rtmusk.f +++ b/src/rtmusk.f @@ -114,13 +114,15 @@ subroutine rtmusk real*8 :: xkm, det, yy, c1, c2, c3, c4, wtrin, p, vol, c, rh real*8 :: topw,msk1,msk2,detmax,detmin,qinday,qoutday real*8 :: volrt, maxrt, adddep, addp, addarea - real*8 :: rttlc1, rttlc2, rtevp1, rtevp2 + real*8 :: rttlc1, rttlc2, rtevp1, rtevp2, sum_rttlc, sum_rtevp jrch = 0 jrch = inum1 qinday = 0; qoutday = 0 det = 24. + sum_rttlc = 0.0 + sum_rtevp = 0.0 !! Water entering reach on day @@ -292,14 +294,15 @@ subroutine rtmusk end if end if rttlc = rttlc1 + rttlc2 + + sum_rttlc = sum_rttlc + rttlc end if !! calculate evaporation rtevp = 0. if (rtwtr > 0.) then - - aaa = evrch * pet_day / 1000. + aaa = evrch * pet_day / (1000.0 * nn) if (rchdep <= ch_d(jrch)) then rtevp = aaa * ch_l2(jrch) * 1000. * topw @@ -338,7 +341,7 @@ subroutine rtmusk end if rtevp = rtevp1 + rtevp2 end if - + sum_rtevp = sum_rtevp + rtevp !! define flow parameters for current iteration flwin(jrch) = 0. flwout(jrch) = 0. @@ -371,6 +374,9 @@ subroutine rtmusk !! volinprev(jrch) = wtrin !! qoutprev(jrch) = rtwtr + rttlc = sum_rttlc + rtevp = sum_rtevp + if (rtwtr < 0.) rtwtr = 0. if (rchstor(jrch) < 0.) rchstor(jrch) = 0. diff --git a/src/rtsed.f b/src/rtsed.f index cf5d4d8..1f8dbb7 100644 --- a/src/rtsed.f +++ b/src/rtsed.f @@ -187,7 +187,7 @@ subroutine rtsed !! Bank erosion rchdy(55,jrch) = 0. !! Channel Degredation - rchdy(56,jrch) = deg2 + rchdy(56,jrch) = deg1 + deg2 !! Channel Deposition rchdy(57,jrch) = dep !! Floodplain Deposition diff --git a/src/smeas.f b/src/smeas.f index 2309669..f7b211f 100644 --- a/src/smeas.f +++ b/src/smeas.f @@ -96,8 +96,10 @@ subroutine smeas end do return -! 5200 format (7x,300f8.3) -! 5300 format (i4,i3,300f8.3) - 5200 format (7x,1800f8.3) - 5300 format (i4,i3,1800f8.3) + + !5200 format (7x,1800f8.3) + !5300 format (i4,i3,1800f8.3) + 5200 format (7x,1900f8.3) !! for pouya + 5300 format (i4,i3,1900f8.3) !! for pouya + end \ No newline at end of file diff --git a/src/sw_init.f b/src/sw_init.f new file mode 100644 index 0000000..da36f81 --- /dev/null +++ b/src/sw_init.f @@ -0,0 +1,46 @@ + subroutine sw_init + + use parm + + integer :: ly + real :: dep_prev + character (len=80) :: titldum + character (len=80) :: header + + open (1234,file='sw_data.in') + open (1235,file='sw_data.out') + + write (1235,100) + 100 format (2x,'HRU',8x,'SOIL LAYER - SOL_ST(mm)',/, 12x,' 1 2 3 4 5 + & 6 7 8 9 10') + + read (1234,*) titldum + read (1234,*) header + + do j = 1, nhru + read (1234,*) k, sw_up10, sw_lo10 + dep_prev = 0. + do ly = 1, sol_nly(j) + if (dep_prev < 100. .and. sol_z(ly,j) < 100.) then + sol_st(ly,j) = sw_up10 * sol_ul(ly,j) + end if + if (dep_prev < 100. .and. sol_z(ly,j) > 100.) then + thick = sol_z(ly,j) - dep_prev + sw_up = (100. - dep_prev) / thick * sol_ul(ly,j) + sw_lo = (sol_z(ly,j) - 100.) / thick * sol_ul(ly,j) + sol_st(ly,j) = sw_up10 * sw_up + sw_lo10 * sw_lo + end if + if (dep_prev > 100. .and. sol_z(ly,j) > 100.) then + sol_st(ly,j) = sw_lo10 * sol_ul(ly,j) + end if + dep_prev = sol_z(ly,j) + end do + write (1235,101) j, (sol_st(ly,j), ly = 1, sol_nly(j)) + end do + +101 format (1x,i4,10f12.3) + close (1234) + close (1235) + + return + end \ No newline at end of file diff --git a/src/transfer.f b/src/transfer.f index 1754a24..f5aa887 100644 --- a/src/transfer.f +++ b/src/transfer.f @@ -123,7 +123,10 @@ subroutine transfer use parm integer :: k, ii - real*8 :: volum, tranmx, ratio + real :: volum, tranmx, ratio,vtot + + nhyd_tr = ih_tran(inum5) + vartran(:,inum3) = 0. !! check beg/end months summer or winter if (mo_transb(inum5) < mo_transe(inum5)) then @@ -134,13 +137,12 @@ subroutine transfer !! compute volume of water in source volum = 0. if (ihout == 2) then - volum = res_vol(inum1) + volum = res_vol(inum1) + varoute(2,nhyd_tr) else volum = rchdy(2,inum1) * 86400. end if if (volum <= 0.) return - nhyd_tr = ih_tran(inum5) !! compute maximum amount of water allowed to be transferred tranmx = 0. @@ -157,19 +159,41 @@ subroutine transfer if (tranmx > 0.) then - !! TRANSFER WATER TO DESTINATION -! select case (inum2) -! case (1) !! TRANSFER WATER TO A CHANNEL -! rchstor(inum3) = rchstor(inum3) + tranmx -! -! case (2) !! TRANSFER WATER TO A RESERVOIR -! res_vol(inum3) = res_vol(inum3) + tranmx -! end select - - !! SUBTRACT AMOUNT TRANSFERED FROM SOURCE + !! Source is a reservoir if (ihout == 2) then - res_vol(inum1) = res_vol(inum1) - tranmx + ratio = 1. - tranmx / volum + ratio1 = 1.- ratio + res_vol(inum1) = res_vol(inum1) * ratio !!|m^3 H2O |water + res_nh3(inum1) = res_nh3(inum1) * ratio !!|kg N |amount of ammonia in reservoir + res_no2(inum1) = res_no2(inum1) * ratio !!|kg N |amount of nitrite in reservoir + res_no3(inum1) = res_no3(inum1) * ratio !!|kg N |amount of nitrate in reservoir + res_orgn(inum1)= res_orgn(inum1)* ratio !!|kg N |amount of organic N in reservoir + res_orgp(inum1)= res_orgp(inum1)* ratio !!|kg P |amount of organic P in reservoir + res_solp(inum1)= res_solp(inum1)* ratio !!|kg P |amount of soluble P in reservior + res_chla(inum1)= res_chla(inum1)* ratio !!|kg chl-a |amount of chlorophyll-a leaving reaservoir + do ii = 2, mvaro + varoute(ii,nhyd_tr) = varoute(ii,nhyd_tr) * ratio + end do + + !!save vartran to add in rchinit/resinit + vartran(2,inum3) = tranmx + vartran(3,inum3) = res_sed(inum1) * tranmx + vartran(4,inum3) = (res_orgn(inum1) + varoute(4,nhyd_tr)) + & / ratio * ratio1 + vartran(5,inum3) = (res_orgp(inum1) + varoute(5,nhyd_tr)) + & / ratio * ratio1 + vartran(6,inum3) = (res_no3(inum1) + varoute(6,nhyd_tr)) + & / ratio * ratio1 + vartran(7,inum3) = (res_solp(inum1) + varoute(7,nhyd_tr)) + & / ratio * ratio1 + + vartran(11,inum3) = lkpst_conc(inum1) * tranmx !mg pesticide + vartran(13,inum3) = (res_chla(inum1) + varoute(13,nhyd_tr)) + & / ratio * ratio1 + + else + !! Source is a reach xx = tranmx ! if (xx > rchstor(inum1)) then ! xx = tranmx - rchstor(inum1) @@ -229,17 +253,17 @@ subroutine transfer rchdy(25,inum1) = rchdy(25,inum1) * ratio rchdy(38,inum1) = rchdy(38,inum1) * ratio rchdy(39,inum1) = rchdy(39,inum1) * ratio - end if - !!subtract from source - do ii = 3, mvaro - varoute(ii,nhyd_tr) = varoute(ii,nhyd_tr) * ratio - end do - !!save vartran to add in rchinit and resinit - vartran(2,inum3) = xx - do ii = 3, mvaro - vartran(ii,inum3) = varoute(ii,nhyd_tr) * ratio1 - end do + !!subratct from source + do ii = 3, mvaro + varoute(ii,nhyd_tr) = varoute(ii,nhyd_tr) * ratio + end do + !!save vartran to add in rchinit + vartran(2,inum3) = varoute(2,nhyd_tr) / ratio * ratio1 + do ii = 3, mvaro + vartran(ii,inum3) = varoute(ii,nhyd_tr) * ratio1 + end do + end if end if diff --git a/src/wetlan.f b/src/wetlan.f index ac198c5..563e2af 100644 --- a/src/wetlan.f +++ b/src/wetlan.f @@ -140,7 +140,7 @@ subroutine wetlan integer :: j, iseas real*8 :: vol, cnv, sed, wetsa, xx, phosk, nitrok, tpco real*8 :: wetsani, wetsili, wetclai, wetsagi, wetlagi - real*8 :: san, sil, cla, sag, lag, inised, finsed,setsed,remsetsed + real*8 :: san, sil, cla, sag, lag, inised, finsed,setsed,remsetsed real*8 :: wetsano, wetsilo, wetclao, wetsago, wetlago real*8 :: qdayi, latqi @@ -175,26 +175,14 @@ subroutine wetlan !! calculate water balance for day wetsa = 0. - wetsa = bw1(j) * wet_vol(j) ** bw2(j) + wetsa = hru_fr(j) * bw1(j) * wet_vol(j) ** bw2(j) wetev = 10. * evwet(j) * pet_day * wetsa wetsep = wet_k(j) * wetsa * 240. wetpcp = subp(j) * wetsa * 10. !! calculate water flowing into wetland from HRU - fr_cur = wet_fr(j) * ((hru_ha(j) - wetsa) / hru_ha(j)) - fr_cur = Max(0.,fr_cur) wetflwi = qday + latq(j) - - if (iwetile(j) == 1) then - wetflwi = wetflwi + qtile - qtile = qtile * (1. - fr_cur) - end if - if (iwetgw(j) == 1) then - wetflwi = wetflwi + gw_q(j) - gw_q(j) = gw_q(j) * (1. - fr_cur) - end if - wetflwi = wetflwi * 10. * (hru_ha(j) * wet_fr(j) - wetsa) qdayi = qday latqi = latq(j) diff --git a/src/wmeas.f b/src/wmeas.f index b69a887..d466516 100644 --- a/src/wmeas.f +++ b/src/wmeas.f @@ -95,8 +95,9 @@ subroutine wmeas end do return -! 5200 format (7x,300f8.3) -! 5300 format (i4,i3,300f8.3) - 5200 format (7x,1800f8.3) - 5300 format (i4,i3,1800f8.3) + + !5200 format (7x,1800f8.3) + !5300 format (i4,i3,1800f8.3) + 5200 format (7x,1900f8.3) !! Pouya + 5300 format (i4,i3,1900f8.3) !! Pouya end \ No newline at end of file diff --git a/src/zero0.f b/src/zero0.f index 28b1980..15a3e1d 100644 --- a/src/zero0.f +++ b/src/zero0.f @@ -505,7 +505,7 @@ subroutine zero0 sol_silt = 0. sol_clay = 0. !! added for Srini in output.mgt nitrogen and phosphorus nutrients per JGA by gsm 9/8/2011 - sol_sumn03 = 0. + sol_sumno3 = 0. sol_sumsolp = 0. strsw = 1. strsw_sum = 0. diff --git a/src/zero_urbn.f b/src/zero_urbn.f index d3d9ca8..5aee35b 100644 --- a/src/zero_urbn.f +++ b/src/zero_urbn.f @@ -18,7 +18,7 @@ subroutine zero_urbn !! subdaily sediment modeling by J.Jeong hhsedy=0. - spl_eros = 0. + eros_spl = 0. rill_mult = 0. eros_expo = 0. sig_g = 0. @@ -81,7 +81,7 @@ subroutine zero_urbn dtp_iyr = 0 dtp_numweir = 0 dtp_numstage = 0 - stp_stagdis = 0 + dtp_stagdis = 0 dtp_reltype = 0 dtp_onoff = 0 dtp_evrsv = 0. @@ -146,7 +146,7 @@ subroutine zero_urbn hrnopcp = 100. ri_sed_cumul = 0. irmmdt = 0. - subdr_kg = 0. + subdr_km = 0. subdr_ickm = 0. num_noirr = 0 ri_qloss = 0