Update to 2012Rev591

Author: crazyzlj

VERSIONS

@@ -1,3 +1,3 @@
 VERSION_MAJOR 2012
-VERSION_MINOR 582
+VERSION_MINOR 591
 VERSION_PATCH

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

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

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

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