subbasin.f

水文模型的原始代码

        !! calculate soil temperature for soil layers
        call solt

        if (pot_vol(j) < 1.e-3 .and. imp_trig(nro(j),nrelease(j),j)==1) &
     &                                                            then
          !! calculate surface runoff if HRU is not impounded or an 
          !! undrained depression
          call surface

          !! compute effective rainfall (amount that percs into soil)
          inflpcp = Max(0.,precipday - surfq(j))
        end if
         
        !! perform soil water routing
        call percmain

        !! compute evapotranspiration
        call etpot
        if (pot_vol(j) < 1.e-6) call etact

        !! compute water table depth using climate drivers
        call wattable

        !! new CN method
        sci(j) = sci(j) +  pet_day * exp(-cncoef * sci(j) / smx(j))     &
     &        - precipday + qday

        !! apply fertilizer-check day and heat units
        do while (iida == ifert(nro(j),nfert(j),j))
          call fert
        end do
        if (igro(j) == 0) then
          if (phubase(j) > phun(nro(j),nfert(j),j)) call fert
        else
          if (phuacc(j) > phun(nro(j),nfert(j),j)) call fert
        end if
 
        !! apply fertilizer/manure in continuous fert operation
        call confert
 
        !! apply pesticide in continuous pest operation
        call conapply

        !! remove biomass from grazing and apply manure
        call graze

        !! compute crop growth
        call crpmd

        !! compute actual ET for day in HRU
        etday = ep_day + es_day + canev

        !! compute nitrogen and phosphorus mineralization 
        call nminrl
        call nitvol
        call pminrl

        !! compute ground water contribution
        call gwmod

        !! apply pesticide
        do while (iida == ipst(nro(j),npest(j),j))
          call apply
        end do
        if (phupst(nro(j),npest(j),j) > 0.) then
          if (igro(j) == 0) then
            if (phubase(j) > phupst(nro(j),npest(j),j)) call apply
          else
            if (phuacc(j) > phupst(nro(j),npest(j),j)) call apply
          end if
        end if

        !! compute pesticide washoff   
        if (precipday >= 2.54) call washp

        !! compute pesticide degradation
        call decay

        !! compute pesticide movement in soil
        call pestlch

        if (surfq(j) > 0. .and. peakr > 1.e-6) then
          if (precipday > 0.) then
            call enrsb(0)
            if (sedyld(j) > 0.) call pesty(0)
            call orgn(0)
            call psed(0)
          end if
        end if

        !! add nitrate in rainfall to soil profile
        call nrain

        !! compute nitrate movement leaching
        call nlch

        !! compute phosphorus movement
        call solp

        !! compute chl-a, CBOD and dissolved oxygen loadings
        call subwq

        !! compute bacteria transport
        call bacteria

        !! compute loadings from urban areas
        if (iurban(j) > 0) call urban

        !! compute undrained depression/impounded area (eg rice) processes
        if (pot_fr(j) > 0.) call pothole
        !! Check date for release/impounding water on rice fields
        if (iida == irelease(nro(j),nrelease(j),j)) then
          nrelease(j) = nrelease(j) + 1
        else if (phuimp(nro(j),nrelease(j),j) > 0.0001) then
          if (igro(j) == 0) then
            if (phubase(j) > phuimp(nro(j),nrelease(j),j)) then
              nrelease(j) = nrelease(j) + 1
            end if
          else 
            if (phuacc(j) > phuimp(nro(j),nrelease(j),j)) then
              nrelease(j) = nrelease(j) + 1
            end if
          end if
        end if

        !! compute sediment loading in lateral flow and add to sedyld
        call latsed

        !! compute nutrient loading in groundwater flow
        call gwnutr
        call gw_no3

        !! lag nutrients and sediment in surface runoff
        call surfstor

        !! lag subsurface flow and nitrate in subsurface flow
        call substor

        !! compute reduction in pollutants due to edge-of-field filter strip
        if (filterw(j) > 0.) then
          call filter
          call buffer
        end if

        !! compute water yield for HRU
        qdr(j) = qday + latq(j) + gw_q(j) + qtile
        if (qdr(j) < 0.) qdr(j) = 0.
        if (qdr(j) > 0.) then
          qdfr = qday / qdr(j)
        else
          qdfr = 0.
        end if

        !! compute wetland processes
        call wetlan

        !! compute pond processes
        call hrupond

        !! perform irrigation operations from shallow aquifer, deep
        !! aquifer and sources outside watershed
        if (irrsc(j) > 2) then
          !! irrigation operation
          if (iir(nro(j),nirr(j),j) > 0) then
            if (iida == iir(nro(j),nirr(j),j)) call irrsub
          else
            if (phuacc(j) > phuirr(nro(j),nirr(j),j)) call irrsub
          endif
          !! auto-irrigation operation
          if (auto_wstr(nro(j),nair(j),j) > 0.) call autoirr
        end if

        !! consumptive water use (ponds, shallow aquifer, deep aquifer)
        call watuse

        !! perform water balance
        call watbal

      endif

      !! perform output summarization
      call sumv

      !! summarize output for multiple HRUs per subbasin
      !! store reach loadings for new fig method
      call virtual

      ihru = ihru + 1
      end do

      return
      end