Merge tag 'rev585_vanilla'

Author: mlt

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))

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-qstage<qaddon) qstage = qin - qaddon
-               
-               !! Flow over the emergency weir
-                if (k==dtp_numstage(sb)) then
-                   qstage = qstage + dtp_cdis(sb,k) * 1.84 * 
-     &             dtp_totwrwid(sb) * (qdepth ** 1.5) * 60. * idt !m3/s
-                end if
-                qout = qout + qstage 
-          
-            !! Partially submerged
-	         else    
-	            watdepact = qdepth - dtp_addon(sb,k)  !! look for add on
-	            if (watdepact<0) watdepact = 0.  !! created for add on
-                if (dtp_weirtype(sb,k)==2) then
-                   dtp_wrwid(sb,k) = dtp_diaweir(sb,k) * 0.667
-                end if
-		         !! Estimate weir/orifice discharge
-		        qstage = dtp_cdis(sb,k) * 1.84 * dtp_wrwid(sb,k) * 
+		        
+		          !! 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	         
+               else
+               !! Partially submerged
+                  watdepact = max(qdepth - dtp_addon(sb,k),0.)
+                  dtp_wrwid(sb,k) = dtp_diaweir(sb,k) * 0.667
+		           
+		          !! weir/orifice discharge
+		          qstage = dtp_cdis(sb,k) * 1.84 * dtp_wrwid(sb,k) * 
+     &		            watdepact ** 1.5 !m3/s
+		          qstage = qstage * dtp_numweir(sb) * 60. * idt !m^3
+               end if               
+	         
+	        else
+	        !! Rectangular weir
+	           watdepact = max(qdepth - dtp_addon(sb,k),0.)
+
+               !! Estimate weir/orifice discharge
+		       qstage = dtp_cdis(sb,k) * 1.84 * dtp_wrwid(sb,k) * 
      &		       watdepact ** 1.5 !m3/s
-		        qstage = qstage * dtp_numweir(sb) * 60. * idt !m^3
-
-	            !Limit total outflow amount less than available water above addon
-	            if(qin-qstage<qaddon) qstage = qin - qaddon
-
-	   	        ! Limit overflow amount to the volume above add-on level
-		        if(qin>qaddon) 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.

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.

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 

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)<sp_pd(i,k)) then
-           write (*,*) " "
-           write (*,*) "Urban BMP Warning!!"
-           write (*,*) "In subbasin ", i
-           write (*,*) "The outlet pipe diameter for sandfilter is"
-           write (*,*) " larger than that for sedimentation basin."
-           write (*,*) "The filter pipe diameter was automatically"
-           write (*,*) "corrected to the pipe size of the "
-           write (*,*) "sedimentation basin."
-           ft_pd(i,k) = sp_pd(i,k)
-         endif
+         if (ft_pd(i,k)<sp_pd(i,k)) ft_pd(i,k) = sp_pd(i,k)
 
          if (ft_dep(i,k)<100) ft_dep(i,k) = 100.
          if (sf_ptp(i,k)>1) sf_ptp(i,k) = 1