pond.f

水文模型的原始代码

      subroutine pond(k)
      
!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine routes water and sediment through ponds
!!    and computes evaporation and seepage from the ponds

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bp1(:)      |none          |1st shape parameter for pond surface area
!!                               |equation
!!    bp2(:)      |none          |2nd shape parameter for the pond surface area
!!                               |equation
!!    chlap(:)    |none          |chlorophyll-a production coefficient for pond
!!    hru_sub(:)  |none          |subbasin in which HRU/reach is located
!!    iflod1(:)   |none          |beginning month of non-flood season
!!    iflod2(:)   |none          |ending month of non-flood season
!!    ipnd1(:)    |none          |beginning month of 2nd "season" of nutrient
!!                               |settling
!!    ipnd2(:)    |none          |ending month of 2nd "season" of nutrient
!!                               |settling
!!    i_mo        |none          |current month of simulation
!!    ndtarg(:)   |none          |number of days required to reach target
!!                               |storage from current pond storage
!!    nsetlp(1,:) |m/day         |nitrogen settling rate for 1st season
!!    nsetlp(2,:) |m/day         |nitrogen settling rate for 2nd season
!!    pet_day     |mm H2O        |potential evapotranspiration on day
!!    pnd_evol(:) |m^3 H2O       |volume of water required to fill pond
!!                               |to the emergency spillway
!!    pnd_fr(:)   |none          |fraction of HRU/subbasin area that drains
!!                               |into ponds
!!    pnd_k(:)    |mm/hr         |hydraulic conductivity through bottom of
!!                               |ponds
!!    pnd_no3(:)  |kg N          |amount of nitrate originating from surface
!!                               |runoff in pond at beginning of day
!!    pnd_no3g(:) |kg N          |amount of nitrate originating from 
!!                               |groundwater in pond at beginning of day
!!    pnd_no3s(:) |kg N          |amount of nitrate originating from lateral
!!                               |flow in pond at beginning of day
!!    pnd_nsed(:) |kg/L          |normal ratio of sediment to water in pond
!!    pnd_orgn(:) |kg N          |amount of organic N originating from 
!!                               |surface runoff in pond at beginning of day
!!    pnd_orgp(:) |kg P          |amount of organic P originating from 
!!                               |surface runoff in pond at beginning of day
!!    pnd_psed(:) |kg P          |amount of mineral P attached to sediment
!!                               |originating from surface runoff in pond at
!!                               |beginning of day
!!    pnd_pvol(:) |m^3 H2O       |volume of water required to fill pond
!!                               |to the principal spillway
!!    pnd_sed(:)  |kg/L          |ratio of sediment to water in pond
!!    pnd_solp(:) |kg P          |amount of soluble P originating from surface
!!                               |runoff in pond at beginning of day
!!    pnd_solpg(:)|kg P          |amount of soluble P originating from
!!                               |groundwater in pond at beginning of day
!!    pnd_vol(:)  |m^3 H2O       |volume of water in pond
!!    pndflwi     |m^3 H2O       |volume of water flowing into pond on day
!!    pndsedin    |metric tons   |sediment entering pond during day
!!    psetlp(1,:) |m/day         |phosphorus settling rate for 1st season
!!    psetlp(2,:) |m/day         |phosphorus settling rate for 2nd season
!!    seccip(:)   |none          |water clarity coefficient for pond
!!    sed_stl(:)  |kg/kg         |fraction of sediment remaining suspended in
!!                               |impoundment after settling for one day
!!    sol_sumfc(:)|mm H2O        |amount of water held in the soil profile
!!                               |at field capacity
!!    sol_sw(:)   |mm H2O        |amount of water stored in the soil profile
!!                               |on any given day
!!    subp(:)     |mm H2O        |precipitation for the day in HRU
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    pnd_chla(:) |kg chl_a      |amount of chlorophyll-a in pond at end of day
!!    pnd_no3(:)  |kg N          |amount of nitrate originating from surface
!!                               |runoff in pond at end of day
!!    pnd_no3g(:) |kg N          |amount of nitrate originating from
!!                               |groundwater in pond at end of day
!!    pnd_no3s(:) |kg N          |amount of nitrate originating from lateral
!!                               |flow in pond at end of day
!!    pnd_orgn(:) |kg N          |amount of organic N originating from
!!                               |surface runoff in pond at end of day
!!    pnd_orgp(:) |kg P          |amount of organic P originating from
!!                               |surface runoff in pond at end of day
!!    pnd_psed(:) |kg P          |amount of mineral P attached to sediment
!!                               |originating from surface runoff in pond at
!!                               |end of day
!!    pnd_seci(:) |m             |secchi-disk depth of pond
!!    pnd_sed(:)  |kg/L          |ratio of sediment to water in pond
!!    pnd_solp(:) |kg P          |amount of soluble P originating from surface
!!                               |runoff in pond at end of day
!!    pnd_solpg(:)|kg P          |amount of soluble P originating from
!!                               |groundwater in pond at end of day
!!    pnd_vol(:)  |m^3 H2O       |volume of water in pond
!!    pndev       |m^3 H2O       |evaporation from pond on day
!!    pndflwo     |m^3 H2O       |volume of water flowing out of pond on day
!!    pndpcp      |m^3 H2O       |precipitation on pond during day
!!    pndsedc     |metric tons   |net change in sediment in pond during day
!!    pndsedo     |metric tons   |sediment leaving pond during day
!!    pndsep      |m^3 H2O       |seepage from pond on day
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    chlaco      |ppb (ug/L)    |concentration of chlorophyll-a in pond
!!    iseas       |none          |nutrient settling rate season
!!    k           |none          |HRU or reach number
!!    nitrok      |none          |fraction of nitrogen in pond removed by
!!                               |settling
!!    phosk       |none          |fraction of phosphorus in pond removed by
!!                               |settling
!!    pndsa       |ha            |surface area of pond on current day
!!    sed         |kg/L          |sediment concentration in pond at beginning of
!!                               |day
!!    targ        |m^3 H2O       |target storage level in pond
!!    tpco        |ppb (ug/L)    |concentration of phosphorus in pond water
!!                               |on day
!!    vol         |m^3 H2O       |volume of water in pond at beginning of day
!!    xx          |none          |variable to hold intermediate calc result
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Min

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      integer, intent (in) :: k
      real :: vol, sed, pndsa, xx, targ, tpco, phosk, nitrok, chlaco
      integer :: iseas


        !! store initial values
        vol = 0.
        sed = 0.
        vol = pnd_vol(k)
        sed = pnd_sed(k)

        !! calculate water balance for day
        pndsa = 0.
        pndsa = 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.

        !! new water volume for day
        pnd_vol(k) = pnd_vol(k) - pndsep - pndev + pndpcp + pndflwi

        if (pnd_vol(k) < 0.001) then
          !! if volume deficit in pond reduce seepage
          !! so that the pond volume is zero
          pndsep = pndsep + pnd_vol(k)
          pnd_vol(k) = 0.
          !! if seepage is less than the volume deficit, take the remainder
          !! from evaporation
          if (pndsep < 0.) then
            pndev = pndev + pndsep
            pndsep = 0.
          end if
          pnd_sed(k) = 0.
          pnd_solp(k) = 0.
          pnd_psed(k) = 0.
          pnd_orgp(k) = 0.
          pnd_solpg(k) = 0.
          pnd_orgn(k) = 0.
          pnd_no3(k) = 0.
          pnd_no3s(k) = 0.
          pnd_no3g(k) = 0.
          pnd_chla(k) = 0.
          pnd_seci(k) = 0.

        else

          !! compute new sediment concentration
          pnd_sed(k) = (sed * vol + pndsedin) / pnd_vol(k)

          !! compute outflow
          if (pnd_evol(k) <= 0.) then
            !! all storage over principle is outflow
            if (pnd_vol(k) <= pnd_pvol(k)) then
              pndflwo = 0.
            else
              pndflwo = pnd_vol(k) - pnd_pvol(k)
            end if
          elseif (pnd_vol(k) > pnd_evol(k)) then
            !! if emergency level is defined, anytime pond volume
            !! exceeds this level, all excess is released
            pndflwo = pnd_vol(k) - pnd_evol(k)
          else
            !! target storage based on flood season and soil water
            xx = 0.
            targ = 0.
            if (iflod2(k) > iflod1(k)) then
              if (i_mo > iflod1(k) .and. i_mo < iflod2(k)) then
                targ = pnd_evol(k)
              else
                xx = Min(sol_sw(k) / sol_sumfc(k),1.)
                targ = pnd_pvol(k) + .5 * (1. - xx) * (pnd_evol(k) -    &
     &                                                      pnd_pvol(k))
              end if
            else
              if (i_mo > iflod1(k) .or. i_mo < iflod2(k)) then
                targ = pnd_evol(k)
              else
                xx = Min(sol_sw(k) / sol_sumfc(k),1.)
                targ = pnd_pvol(k) + .5 * (1. - xx) * (pnd_evol(k) -    &
     &                                                      pnd_pvol(k))
              end if
            end if
            if (pnd_vol(k) > targ) then
              pndflwo = (pnd_vol(k) - targ) / ndtarg(k)
            else
              pndflwo = 0.
            end if
          end if

          !! compute final pond volume
          pnd_vol(k) = pnd_vol(k) - pndflwo
          if (pnd_vol(k) < 0.) then
            pndflwo = pndflwo + pnd_vol(k)
            pnd_vol(k) = 0.
          endif

          !! compute change in sediment concentration due to settling
          if (pnd_sed(k) > pnd_nsed(k)) then
            pnd_sed(k) = (pnd_sed(k) - pnd_nsed(k)) * sed_stl(k) +      &
     &                                                       pnd_nsed(k)
          end if

          !! compute sediment leaving pond
          pndsedo = pnd_sed(k) * pndflwo
 
          !! net change in amount of sediment in pond for day
          pndsedc = vol * sed + pndsedin - pndsedo - pnd_sed(k) *       &
     &                                                        pnd_vol(k)

          !! determine settling rate
          !! part of equation 29.1.3 in SWAT manual
          if (i_mo >= ipnd1(k) .and. i_mo <= ipnd2(k)) then
            iseas = 1
          else
            iseas = 2
          endif
          phosk = 0.
          nitrok = 0.
          phosk = psetlp(iseas,k) * pndsa * 10000. / pnd_vol(k)  !setl/mean depth
          phosk = Min(phosk, 1.)
          nitrok = nsetlp(iseas,k) * pndsa * 10000. / pnd_vol(k) !setl/mean depth
          nitrok = Min(nitrok, 1.)

          !! remove nutrients by settling
          !! other part of equation 29.1.3 in SWAT manual
          pnd_solp(k) = pnd_solp(k) * (1. - phosk)
          pnd_psed(k) = pnd_psed(k) * (1. - phosk)
          pnd_orgp(k) = pnd_orgp(k) * (1. - phosk)
          pnd_solpg(k) = pnd_solpg(k) * (1. - phosk)
          pnd_orgn(k) = pnd_orgn(k) * (1. - nitrok)
          pnd_no3(k) = pnd_no3(k) * (1. - nitrok)
          pnd_no3s(k) = pnd_no3s(k) * (1. - nitrok)
          pnd_no3g(k) = pnd_no3g(k) * (1. - nitrok)

          tpco = 0.
          if (pnd_vol(k) + pndflwo > 0.1) then
          tpco = 1.e+6 * (pnd_solp(k) + pnd_orgp(k) + pnd_psed(k) +     &
     &                            pnd_solpg(k)) / (pnd_vol(k) + pndflwo)
          else
            tpco = 0.
          endif
          chlaco = 0.
          pnd_chla(k) = 0.
          pnd_seci(k) = 0.
          if (tpco > 1.e-4) then
            !! equation 29.1.6 in SWAT manual
            chlaco = chlap(k) * 0.551 * (tpco**0.76)
            pnd_chla(k) = chlaco * (pnd_vol(k) + pndflwo) * 1.e-6
          endif
          if (chlaco > 1.e-4) then
            !! equation 29.1.8 in SWAT manual
            pnd_seci(k) = seccip(k) * 6.35 * (chlaco**(-0.473))
          endif
        end if

      return
      end